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1. Introduction. Although the propagation of sound waves in moving media has 
received considerable attention [1, - - - , 9],! little information is available concern- 
ing the propagation of such disturbances in rotational streams or concerning the 
propagation of transient rotational phenomena. It is shown in the present paper that 
the wave fronts associated with those parts of a disturbance which are derivable from 
a potential propagate in a rotational stream according to those laws which they are 
already known to obey in an irrotational stream.’ It is further shown that the rota- 
tional disturbances drift with the stream rather than propagate relative to the moving 
fluid. 

The analysis consists of an application of conventicnal perturbation procedures 
to the Navier-Stokes and continuity equations. The equations so derived are treated 
according to the theory of characteristics. The results obtained lead to a general ex- 
pression for the Mach lines of an arbitrary supersonic flow and also suggest a new 
method of wind tunnel calibration which eliminates the necessity of placing an ob- 
stacle in that portion of the stream being calibrated. Finally, predictions are carried 
out as to the nature of pulses which are formed at a surface and then propagate 
through a boundary layer into a uniform stream. 

2. The equations of motion. In this analysis, we shall consider the propagation 
of small disturbances in fluid streams which are characterized by three functions of 
the space coordinates and the time, namely: po (the density), po (the pressure), and 
Vo (the velocity). No restrictions will be applied to these functions except that they 
obey the differential equations implying the conservation of momentum, mass, and 
energy. These equations, known familiarly as the Navier-Stokes, continuity, and en- 
ergy equations, may be written in the forms: 


1 
(v-grad)v + ov/at + — grad p = — L(v) (1) 
p p 


div v + @ In p/dt + v-grad In p = 0. 


dU/dt + pd(p)/dt = Q + x 


* Received Jan. 10, 1946. 
1 Numbers in brackets refer to the bibliography. 
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In the foregoing equations, y is the viscosity of the fluid, L symbolizes (A+ grad div) 
where A is the Laplacian operator; Q is the rate of heat accumulation, U is the internal 
energy, and x abbreviates the viscous dissipation terms. Discussions of these equa- 
tions are conveniently found in [7] and [8]. The necessity of manipulating the energy 
equation in the investigation may be eliminated by using the following assumption. 
The changes in pressure and density accompanying the disturbance are taken to obey 


the law 

b/ po = (p/p) (3’) 
where p, ~, Vv, characterize the disturbed stream; that is, the disturbance is a phe- 
nomenon such that the changes in state from undisturbed to disturbed stream are 
isentropic. Note that this in no way restricts fo and po. The appendix indicates briefly 
the fact that while this assumption is by no means rigorously justified, it leads to valid 
results. 

It is convenient at this point to introduce the small parameter e. Although this 
may be done in a fairly arbitrary manner, we shall define it in the following way in 
order to avoid any possible ambiguities. Let the initial conditions of any particular 
problem be such that at time zero, p = po+€p1, where the maximum value of p1/po over 
the region under consideration is unity. Thus, since we are considering small disturb- 
ances, € is a small number compared to unity. Consistent with this notion, we shall 
write p=potepit +e7pet - ++, P=potepfit ---,and v=vo+evit ---, at time 
t; and shall require that the series be valid over a range of ¢ Since disturbances can 
usually be expected to attenuate, it is certainly reasonable to expect that the series 
will converge for sufficiently small values of this parameter. 

If we now substitute the foregoing forms of p, p, and v, into Eqs. (1) and (2), 
eliminate the p; (except for po) by using Eg. (3’), and collect the coefficients of each 
power of e€, we obtain: 


B 


1 
(vo-grad)Vo + dvo/dt + — grad po L(vo) 
Po 


Po 


+e {(wo-srad)y; + (v,-grad)vo + dv,/dt + grad (<: =) 


Po 
p 2?p Bb p 
mS ont f+ he ged oe ~ “| x0) om Liv) } a a 
p* Po Po Po f 
and 
0 In po/ dt + div Vo + Vo grad In Po 
‘ Pl 7) Pi a 
+ €4div v; + vi: grad In po + Vo: grad ——- + —( — +---=0Q, (5) 
Po Ot Po 


where, a0 =Yfo/po.- 
When, in equations (4) and (5), the coefficients of e® are equated to zero, we find 


two of the necessary conditions that the functions po, po, Vo, characterize a possible 
fluid stream. Since these quantities must vanish identically, we may omit them from 
Eqs. (4) and (5), and divide the remaining equalities through by e. When we allow 
€ to approach zero, we see that the mathematically exact solution to the problem is 
found by setting the coefficients of € in Eqs. (4) and (5) to zero, Hence, we may expect 
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that the functions pi, vi, so determined will provide a good first approximation to the 
behaviour of small amplitude disturbances. This, of course, is the conventional per- 
turbation reasoning. 

If we had been willing to assume at the outset a functional relationship p= p(p) 
applicable both to the stream and the disturbance, the perturbation procedure would 
have been unnecessary. The forthcoming techniques could have been applied directly 
to Eqs. (1) and (2). However, the solution possesses the desired generality only when 
we refrain from such restrictions on the nature of the stream. This leads to a choice 
between working with the energy equation or using the foregoing procedure; the lat- 
ter seems more convenient. As a matter of fact, some of the results of this analysis 
differ from those of previous investigators only in that they are obtained for any 
stream wherein the medium behaves as a continuum rather than one of a very re- 
stricted character. 

Recalling now that any vector may be expressed as the gradient of a scalar plus 
the curl of a vector and that 


(B-grad)C + (C-grad)B = grad (B-C) + (curl B X C) + (curl C X B), 


one may write 


v,; = grad ¢ + curl A, 


and the differential equations defining p; and vi become 


2p dp ) p 
grad | (wo-¥) + do bed aa ~] a re curl A + as = grad In po 


Po at Po 
p rm p 
— — grad po + @: X Vo + wo X Vi = “| x) + * rive) | (6) 
p> Po Po 
and 
d {pi 
Ag + — (“) + vi-grad In po = 0 (7) 
dt \ po 


where w;=curl v;, and d/dt = [vo- grad +0/0t]. When the operation “curl” is performed 
on Eq. (6), the following equality arises: 


o “a p 2p 
— = curl j is | 0) + * Low) = — grad In po 
ol Po Po Po 
p 
+ = grad po — @1 X Vo — @o X nt ‘ (8) 
a 


It is evident by inspection of Eq. (8) that an identically vanishing initial choice of w; 
does not imply that this function will vanish for all time, as for example, is the case 
in an irrotational stream. Thus, we cannot omit w: in this investigation. 

It will prove useful to define an (artificial) auxiliary potential y in the following 
manner? (Eq. (6) implies the existance of this quantity). 


2 This will allow us to eliminate p;/po and thus obtain equations in which each unknown has the di- 


mensions of a velocity potential. 
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oy OA ok p 
— grad — + curl — = “| 103) + a Liv) | — — grad In po 
ot ot Po Po Po 


p 
+ — grad Po — @1 X Vo — Oo X Vi = Oz (9) 
p? 


0 
Upon substitution of Eq. (9) into Eq. (6), the latter becomes 


2p: Ob @ 
grad | vow + @g— + — + ~ | = 0 
Po Ot ot 
This however, we may solve for p:/po arbitrarily choosing the “constant of integra- 
tion” to be zero.’ We obtain 


2 d 
P1/po = — do | ve-cur A + ° ~|. (10) 
dt ot 
This may be combined with Eq. (7) to give 
d/1 de dt 1 /ay 
Ad — —(— —])+ vi -grad In = =| (= + ve-curt a) | (11 
' < 4 ia ee dt La Nae" 


It will be shown directly that the wave front propagation can be derived from 
Eqs. (8), (9), and (11), provided we can justify the omission of the term (u/po)L(v:) 
from Eqs. (8) and (9). We note, considering Eq. (8), that if only the terms 0w;/0/ 
and curl (u/po)L(v1) were non-vanishing, we would have virtually the equation for 
the conduction of heat, that is 0w;/0t=(u/po)Aw:. The “conduction coefficient” is very 
small (in air the spreading of vorticity is known to be very slow) so that the term in 
question may be thought of as one which causes a small dispersive effect. It is to be 
understood, then, that this effect is to be superimposed on any results which are ob- 
tained by treating the equations from which this term has been omitted. With this 
omission we are now ready to apply the method of characteristics. 

Hadamard [1] has shown the following facts concerning second order differential 
equations which will be useful in the analysis of the foregoing equations. He con- 


siders the equation 


n 


D aupi +h = 0 

i, keel 
in the m independent variables x, - - - , Xx», where the a;, and h are functions of the 
unknown quantity 2, the x;, and the first partial derivatives of z with regard to the x;; 
pix =0°2/0x,0x,.. The differential equation which defines the characteristic surfaces 
(wave fronts) of this equation is given by 


n—1 


n—l 
B= > auPsPi — D> ainPi + dan = 0 (12) 


t,k=1 t=1 
where P;=0x,/0x; when the “surface” is written in the form 
x, = Xn(%1, es Sunt). (12a) 


8 Any such (actually time dependent) constant could be absorbed in d¥/d¢ and would contribute 


nothing to ¥. 
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Furthermore, let there be s unknown functions 2, - - - , 2,, and s equations of the form 
Do ainpin + bingix + +++ + Cingie + h = 0. 
i,k 
Here the a,x, biz, + - + are respectively the coefficients of the second derivatives of 
2, 2, ° - -. Lhe characteristic surfaces of this system of equations are determined by 


the relation 


| Bu, Biz, - ++, Br | 


| Bay 
i = 0. (13) 
| Bu Bn 
The B,, are analogous to the quantity B of Eq. (12). In fact, when a takes the values 
1, 2,---+, s, B is derived respectively from the first, second, - - - sth, equations. 


When 8 takes these values, the B,g are obtained from the aix, Diz, - - + , Cix, Tespec- 
tively. The present problem deals with the five unknown quantities ¢, y, and the three 
components of A. Eq. (9) is equivalent to four scalar equations‘ if we specify (for 
example) that y and A are to be those solutions for which div A=0. This is no restric- 
tion since only curl A appears in v;. We may, then, apply the foregoing type of analy- 
sis to Eqs. (9) and (11) (with po/pi replaced by the expression given in Eq. (10)). 
In fact, in order to determine the characteristic surfaces which define a motion involv- 
ing the function ¢, we need only a brief inspection of Eq. (9). It is evident that no sec- 
ond derivatives of @ appear in this equation. Thus, when it is split into its four 
subdivisions, we find that the four quantities Ba, - - - , Bs:, which appear in the left 
column of Eq. (13), vanish. This implies that Eq. (13) is satisfied when either By or 
the minor associated with this quantity vanishes. Since the vanishing of the former 
involves only the coefficients of derivatives of ¢, we may assume that this surface 
will be associated with the potential type of disturbance. The vanishing of the minor 
will correspond to the propagation of disturbances of the rotational type. 

If we now compute By using, of course, the a;, of Eq. (11), we find the same wave 
front equation which was found by Hadamard for the isentropic stream. That is, 
the time-position correlation of a wave front does not depend on the character of the 
stream but only (as the following equation will show) on the local values of the 
quantities %, Yo, Wo, and do. The first three of these are the components of vo. This 
wave front equation, in a form somewhat more convenient for our purposes than 


Hadamard’s, is shown below. 
dy/dt + upAdy/Ax + wody/dz — vo + ao[1 + (Ay/dx)? + (Ay/dz)2]/2? = 0. (14) 


Eq. (8) indicates that whenever w» and ¢ are each non-vanishing in a given region, 
a rotational motion @, is generated continuously. This being so, there is always a 
possible “vorticity wave front” coincident with the wave front associated with @¢. 
Hence, if we treat Eq. (8) according to the foregoing method, using the components 
of curl A as the unknown function, we find that the determinant vanishes identically. 
This is to be expected since the operation which led to Eq. (8) eliminated the higher 
derivatives of @ while retaining the higher derivatives of A. Hence, formally, the char- 


4 Any equation of the form curl M+grad Q=C can be reduced to the forms grad Q=P and curl 
M =N if one is ingenious enough to separate C into the required parts. 
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acteristics method fails to give the desired information. We note that in this‘method, 
however, the only terms which affect the positional nature of the propagation are 
those containing second derivatives of the unknown functions. If, in Eq. (8), we 
segragate these terms of the required order we find that they comprise exactly the 
single term dw;/dt. Therefore, in so far as the position of the disturbance is concerned, 
we have dw; /dt =0; that is, the time rate of change of vorticity, relative to an observer 
moving with the particle, vanishes. In other words, the rotational disturbance drifts 
with the stream instead of propagating relative to it. This statement must be modi- 
fied, of course, by the results of the diffusion-implying viscous term which was omitted 
in this analysis. 

3. The two-dimensional problem. Since, in general, the functions po, po, Vo asso- 
ciated with any given stream are not known (even approximately in many cases), it 
seems of interest to describe a method of wind tunnel calibration based on the fore- 
going analysis (in particular on Eq. 14). This proposed procedure will be seen to have 
the advantage that it does not require the insertion of an obstacle into that portion 
of the stream being calibrated. Let us consider only tunnels which are bounded by 
the side walls z= +b, where b is some constant. In this two-dimensional wind tunnel, 
the flow in the neighborhood of z=0 is essentially independent of z. Let us also restrict 
our consideration to disturbances having reflective symmetry about the plane z=0. 
Then at z=0, Eq. (14) reduces to 


dy/dt + ao[1 + (Ay/dx)?]}/? — v9 + updy/dx = 0. (15) 


We now have an equation, linear in the three quantities which we wish to deter- 
mine; Uo, Yo, and do. Suppose we generate pulses at several points along some bound- 
ary of the stream, say by the use of an electric spark. The wave fronts of these pulses 
may be observed (photographed) at successive time intervals. The values of dy/dx 
and 0y/dt can be determined from the photographs for each pulse throughout the 
region it traverses. For each point traversed by at least three pulses, we may form 
three simultaneous equations in the unknown quantities by using these experimen- 
tally determined values as coefficients in Eq. (15). Figures 1 to 4 illustrate such photo- 
graphs of sound pulses in a fairly uniform stream of air. The development of the 
techniques used in obtaining these Schlieren photographs should be credited to the 
authors of [9]. In [9] the details of the experimental procedure are explained quite 
fully. 
For an isentropic region of the stream (where the stagnation condition is known) 
only two pulses are needed since (see [11 ]) 


ao = (as) -y=- 1)(uo + 0) /2. 


Finally, for an essentially one-dimensional stream (e.g., a jet or a slowly converg- 
ing channel) Eq. (15) becomes 
[2¢y + 1) — 2y — tut]? — 2y 
y+1 


(15a) 





V0/dst = 


where p =a, 'dy/dt. 
When the stream is supersonic, we may generate stationary disturbances (Mach 
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lines) by placing very small irregularities in the boundaries of the passage. For this 


case Eq. (15) reduces to 
2 2 1/2 
— U% taiM —1 
ay/ax = —— ol (16) 
a? — u2 
0 0 





where M, the Mach number, is given by M?=(u2+1)/a3. A mesh of Mach lines pro- 
ceeding from both edges of the passage give sufficient information to calibrate any 
stream known to be isentropic with known stagnation condition. For a one-dimen- 
sional stream, we have the familiar formula for the Mach angle @ 


dy/dx = tan 0 = (M? — 1)", (17) 


Note that when one wishes to find the characteristics of a given supersonic stream, 
the classical Charpit procedure will always provide solutions for Eq. (16). The equa- 
tion analogous to Eq. (16) in three dimensions follows directly from Eq. (14) by 
merely dropping the time dependent term. 

4. The effect of boundary layers. A problem of considerable interest arises in con- 
nection with the ideas of the foregoing section when we inquire into the effect of the 
boundary layer on the form of the wave front when the pulse is generated at the 
surface of a boundary (or obstacle). We shall use the Charpit procedure to solve 
Eq. (15) for this case, using, of course, an idealized group of values for uo, vo, and apo. 
The justification of the steps of this procedure are given in [10] and need not be given 
here, so we shall proceed formally with this method. The solution obtained will be in 
closed form and is readily verified.(as a solution of Eq. 15) by mere substitution. 

We characterize the stream by the functions 


uy =0, vm = vx/6, for «6, vm =v for x26, a = a = const. 


This simplification of an actual situation is somewhat drastic but useful information 
results. We first set £=x/6, n=y/6, r=at/5, M=v/a, p=0n/0E, g=0n/0r, and Eq. 
(15) in the notation of [10] becomes 


Fié,n7290= qt [1+ 97]? - ME=0 for OSES 1 


or 








F=q+ [1+ p?]?—-M=0 for 1<¢. (18) 
We proceed by considering the associated ordinary differential equations 
1 d d d dr 
ap - q i 7 _ g Ps (19) 
F;+ pF, F, + gF, — pF, — gFy —F, —F 


and choose any solution which expresses p or g in terms of a parameter a. In our case, 


(formally) 


a ( .. +0) y= (19a) 
M Tr q/ “7 [1 + p2]1/2 q 7 a 


and g=a is the required solution. When this is substituted into Eq. (18), we obtain 
for p in the respective regions 


pb = [(a — Mt)? —1]'/? and p = [(a — M)? — 1]"/2, (20) 
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We now determine 7 by the following integration 
7—-B= foae + qdr (21) 


8 is chosen to suit the initial and boundary conditions; a is determined by the rela- 
tion 0n/da=0, and the sign of » must be taken consistent with that portion of the 


wave front under consideration. 
For the initial condition 7 =0 at £=7=0 (a point source), and for <1, we have 
for that downstream portion of the wave for which p30, 





n = ar + Bla — ME) — Bla) (22) 
where 
B(a) = ala? — 1]"/? — arc cosh a (23) 
7 Mt Mr 4 1/2 
ar a ae ee 


This solution is valid when £5 S&S min (r, 1); 
fo = [72 ++ M-*]!/2 — 1/M. 


When ££, p20, and we must find 7 by writing 


gE 
n(r, £) = (7, £0) + pdé 
§> 
which results in the formula 


n(r, £) = ar — Bla) — B(a — Mé) (25) 
where @ and 8 are still the quantities defined by Eqs. (23) and (24). 
When we consider the upstream portion of the wave front, a must be negative. 
By the foregoing procedure we obtain 











a oo = B(— ai) + B(ME — a1) (26) 
Mt Mr 4 ass 
2 2 7? — ¢2 
or ; 
mi(&, T) aia n(é, T) + Mér. (28) 


This equation is again valid only when £2 &p. 

We now consider that portion of the wave exterior to the boundary layer (i.e. 
£>1). We must again extend the integration of Eq. (21), this time into the uniform 
stream. Using now p= — [(a—M)?]"/?, we obtain for the downstream portion 


e. 
n= art Bla — Mf) — Ala) — f [a — Mt — 1]¥"4¢ 
= ar + B(a — M) — B(a) —[ (a — M)? — 1]/2(E — 1). (29) 


When 07/da is equated to zero we find 


[(@ — M)? — 1]? 1 1 \ 
—1i= r+ — [(a — M)? — 1]!/2 — — [a? —1]"/2. (30 
— a—M y+zl = 0 i! (30) 
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Eqs. (29) and (30) may be considered as parametric equations for — and y in terms of 
the parameter a for each value of r. In a similar manner we obtain for the upstream 





portion 
m = a7 — B(— a) + B(M — a) + [(a1 — M)? — 1]"/7(E1 — 1) (31) 
and 
[(m — a)? _ 1}1/2 1 2 1/2 1 2 1/2 
-—-i= — —<7r+—la,-—1 — — |(M —a;) -1 ; 32 
bs rt led — 1 2 far - ay" — 1". 30 


In Eqs. (29) to (32), 1+Msa<o and —12a,>-—o. 
The solution is now complete except for surface reflections. Note that the velocity 
at which the point of contact between boundary and wave front moves is 





an! M272 J1/2 
nee = a(t, 0) = E + | 
OT t=_0 
or 
oy vf? /2 
a elttagl - 


Figure 5 illustrates the results of the foregoing section for a stream with Mach 
number .50. The peculiar behavior of the solution for 7 > 4/5 leads us to investigate 
the rays of the propagation. It has been shown [1, 2, 3] that Eq. (14) implies that the 
rays be defined by dx = (lap +uo)dt, dy = (mao+10)dt, and dz = (nao+wo)dt, where |, m, 
n, are the direction cosines of the outwardly directed wave front normal with the 
coordinates axes. In particular, it has been shown that for the conditions prevailing 
in the boundary layer specified above, the rays are given by (see [9]) 


1 = — >> [are cosh (mo — ME) + (mo + M2){(mo — MB)? — 1}"2] (34) 


where mp is the value of m~! at the origin. The value of mo for which the ray becomes 
tangent to the line £=1 is given by mp =1+M. However, any ray associated with a 
uniform stream which is directed parallel to that stream will maintain this orienta- 
tion. Hence, this ray which just becomes tangent to the uniform stream bifurcates 
into the curves shown in Fig. 5. Note that no ray which is once reflected back into the 
boundary layer will ever leave this region. This implies that a sizeable portion of the 
energy of such pulses never leaves the boundary layer. Furthermore, as one can 
readily see from the few rays plotted in the figure, the reflections occur in such a 
manner that interference as well as the extremely turbulent conditions in such a 
region make the observation of waves in such regions improbable. This is borne out 
in Figs. 1 to 4. One large discrepancy between these pictures and the theory is easily 
noticed. The limiting upstream ray given by the theory does not agree too well with 
the evidence of Fig. 4. This is due to the fact that the pulse used in the experiment 
started as one of finite amplitude (probably traveling initially at a speed of about 2a»). 
This means that at first the wave travels as though it were a small amplitude wave in 
a slower stream; hence, less distortion from the shape which would be expected with 
no boundary layer (a family of semi-circles) is the logical result to expect. This, of 
course, is consistent with the observations. One other remark is essential in view of 
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the initial idealization of the boundary layer. In the actually occurring physical situa- 
tion, the boundary layer thickens in the direction of flow and thus the sharp break 
in wave front predicted in this theory is not valid. However, the transition from large 
velocity gradient to uniform stream occurs in a sufficiently small region so that little 
energy transfer from the stream part of the wave to the boundary layer is to be ex- 


pected. 





Fic. 5. Predicted wave fronts for the stream defined in section (4). Wave fronts are illustrated for r=.1 
2, /5, 3, 5. The dotted curves are rays of the propagation. The sound source is at the origin. 


If one wishes to account for the large pulse velocity in a mathematical manner, 
he can replace the constant value of a» by a function of the time, large at time zero but 
rapidly approaching the steady value. This, of course, makes the calculations tedious. 

5. Intensity distribution. In general, it is difficult to obtain a solution for ¢ from 
Eq. (11). However, one very interesting solution for the isentropic uniform stream has 


recently appeared. Rott [6] has shown that when 1 =wo=0, uo =const = — Mao, Eq. 
(11) has a solution of the form 
c Mx+R 
¢1 = — cos o( _ a). (35) 
R ao(1 — M?) 


where R= [x?+(1—M7?)(y?+2?) ]!/*. This solution implies a continuous point source 
at the origin and, of course, assumes no boundary layer if the plane (say) y=0 is to 


be a boundary. 
The surfaces of constant phase (characteristic surfaces) are given by 


$ — —————- = const. = A (36) 


a relationship which can be shown to satisfy Eq. (14) for the given stream. 
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In connection with our pulse problem, we see immediately that we may form a 
new solution as the integral 


@ C(w) Mx+R 
= w {| ¢ — ——————_ : 37 
¢ J > cos ( se =) dw (37) 





If Cis properly chosen, this solution corresponds to a pulse of any desired wave form 
originating at the origin. The surfaces of constant phase are again given by Eq. (36) 
and, for a given value of ¢, are circles with centers at the points x= Maot, y=z=0, 
and radii dof. 

Suppose we now choose two values of \ (say A; and dg) to represent two character- 
istic surfaces whose separation (along a radial line from the origin) we can call the 
wave length of the pulse. Then that upstream portion of the pulse at y=z=0 will 
have a wave length (Ae—Ai)ao(1—M) and the downstream section at y=z=0 will 
have a wave length (A,—A1)a0(1+M). That is, the ratio of the thicknesses of these 
two extreme portions of the pulse is (1+M)/(1—™). 

Inspection of Eq. (35) also leads to the conclusion that at these portions of the 
pulse the amplitudes vary in the ratio (1—_M)/(1+M). Thus the amplitude gradients 
at these two sections are in the ratio (1— M)?/(1+™)?. Since the density gradient is 
essentially the quantity observed in the Schlieren optical system, the foregoing con- 
stitutes an explanation of the far superior clarity of the wave front definition in the 
upstream portions of the photographs. This argument has assumed that the pulse 
started at time zero and has a small thickness to radius ratio at the time of observa- 
tion. 

Specifically, the amplitude ratio (i.e., the ratio of amplitude at any point on the 
wave front divided into the amplitude at the upstream extremum) is given by 


Mx 


ates ve 8 
aot(1 — M) - 


Amp. ratio = 1+ M — 


APPENDIX 


We wish to justify here the use of Eq. (3’). We write the energy equation in the 


form 
R aT dp! a i 
Sopa tape Sy (3) 
y-—1 dt dt p p 





where T is the temperature (p=pRT) and x involves a sum of products of the form 
(du/dy)(dv/Ax), (@u/dy)*, - - -. When, as in the foregoing work, the perturbation 
procedure is applied, the terms with coefficient € can be written 


1 d(T/T q 
Soe. Sa ee See (3a) 
: ee 1 dt dt poRT Po 








where we have grouped the terms not containing derivatives of the unknown func- 
tions pi, fi, Vi, 71 in K. Such terms can obviously contribute nothing when the char- 
acteristics method is applied. In Eq. (3a) (y—1)a/R is a ratio of specific heat to 
thermal conductivity, and x: involves products of the form (0u/dy)(00:/0x), 


(Ouo/dy)(Ou1/Oy),- °°. 
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If we can show the third and fourth terms of this equation to be negligible, in- 
tegration of Eq. (3a) leads to the terms with coefficient ¢€ in the expansion of Eq. (3’). 
Using a dimensional treatment analogous to Prandtl’s boundary layer analysis, we 
define a typical length / for the disturbance (say the wave length, if a continuous 
wave is considered, or the breadth of a pulse, etc.) and compare first the terms 
(y—1)-!vo-grad(7i/To) and (u/po)(Ouo/dy)(Ou1/0y). We note that (y—1)~'vo-grad 
(T1/To)~(y—1)-(T1/To)(| v0| /2) and 


MB OU OU; rm | vo| | vi} 
Po oy oy aopo 6 l 


so we have the requirement, if the latter term is to be omitted, that podao/p| v1 | 
>T)/Ti. Since 71/7» and v1 /do are of appreciable magnitude (of order unity) in 
the same region, this inequality is essentially podao/u>>1. Here, 6 is the boundary layer 
thickness or other typical dimension of the stream. Thus we see that except for very 
restricted regions the above inequality will hold and the term in question may be 
omitted. The other terms in x: admit a similar treatment. 

We now compare the terms 0(71/T»)/dt and aAT;/RT po. Note that these are 
exactly the terms which would need to be compared in the still air case; that is, no 
functions characterizing the stream appear except the slowly varying aj. Using the 
typical length / and the fact that the propagation occurs at essentially velocity ado, we 
obtain aAT,/RTopo~aT;/ RT ol?po<KO(T1/T0)/0t~aoT /1T9 as a necessary condition. 
That is, the number a/aolRpo1. This, of course, is the actual situation* and hence 
both terms in question may be omitted. This number is essentially the reciprocal of 
a Prandtl number-Reynolds number product. 
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THIN CYLINDRICAL SHELLS SUBJECTED TO 
CONCENTRATED LOADS* 


BY 


SHAO WEN YUANT 
California Institute of Technology 


Abstract. A single differential equation of the eighth order in the radial displacement is given for 
the equilibrium of an element of a cylindrical shell undergoing small displacements due to a laterally 
distributed external load. The radial deflection of thin cylindrical shells subjected to concentrated, equal 
and opposite forces, acting at the ends of a vertical diameter, is analyzed by the Fourier method. Applica- 
tions of the solution of the problem of the infinitely long cylinder to the problems of a couple acting on 
an infinitely long cylinder in the direction of either the generatrix or the circumference are also discussed. 


1. Introduction. The bending problem of an infinitely long cylinder loaded with 
concentrated, equal and opposite forces, acting at the ends of a vertical diameter, is 
discussed first. The equations of equilibrium of an element of a cylindrical shell un- 
dergoing small displacements due to a laterally distributed external load are reduced 
to a single differential equation of the eighth order in the radial displacement. In 
this equation the various terms are compared as to the order of magnitude and it is 
found that some of the terms are negligible. 

The specified loading function is represented by a Fourier integral in the longitu- 
dinal direction, and by a Fourier series in the circumferential direction. The integral 
representation has the advantage that the boundary conditions are automatically 
taken care of, and no subsequent determination of Fourier coefficients is necessary. 
The Fourier coefficients and the undetermined function in the Fourier integral in this 
case are determined simply from the loading condition. The radial displacement is 
represented in a like manner with the aid of an undetermined function which is ob- 
tained by substituting both radial displacement and loading expressions in the differ- 
ential equation. The definite integrals involved in the expression for radial deflection 
are evaluated by means of Cauchy’s theorem of residues. 

The problem of the inextensional deformation of cylindrical and spherical shells 
was treated in detail by Lord Rayleigh in his “Theory of sound.” The assumption 
of this type of deformation underlies the solution of many problems of practical im- 
portance, such as the determination of stresses in thin cylindrical shells subjected to 
two equal and opposite forces acting at the ends of a diameter or to internal hydro- 
static pressure. It is found that the results obtained in the case of inextensional de- 
formations correspond only to a first approximation of the solution in this paper, and 
the stresses in the proximity of the points of application of the forces are not given 
with sufficient accuracy. 

The expression for the radial deflection of a thin cylinder of finite length is ob- 
tained from the corresponding solution for an infinitely long cylinder by using the 
method of images. It is seen that the difference of these two radial deflections can be 
given by a correction factor included in the expression for a cylinder of finite length. 

* Received April 20, 1945. 

Tt Now at Polytechnic Institute of Brooklyn. 
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The difference is believed to result from restraining the edges at the two ends of the 
finite cylinder. The results indicate that the radial deflection of an infinitely long cyl- 
inder has a very long wave length along the generatrix; however, the wave length 
decreases as the ratio of radius over thickness decreases. It is believed that the long 
wave length phenomenon is due to the elastic reaction along the circumference of the 
shell which can be explained by the radial deflection along the circumference. 
Deflection curves of cylindrical shells with various lengths are calculated and the 
results show that the maximum radial deflection occurs at length over radius ratio 
1/a~™~20. The radial deflection of an infinitely long cylinder with the radius over thick- 
ness ratio a/h=100, becomes zero at about x/a=15 and then reverses its sign. The 
edges of the corresponding cylinder with finite length are so restrained that the nega- 
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Fic. 1. Forces and moments on element of wall. 


tive deflection portion of the infinite cylinder is brought to zero at the edges of the 
cylinder with finite length. Hence, the maximum deflection of a cylinder with //a~20 
is greater than that of the corresponding infinitely long cylinder. 

The problems of a couple acting on an infinitely long cylinder in the direction of 
either the generatrix or the circumference are also analyzed by using the correspond- 
ing solution for the radial deflection under a concentrated load. The action of the 
couple is equivalent to that of two equal and opposite forces acting at an infinitely 
small distance apart. 

2. Fundamental equations. The fundamental equations of a cylindrical shell un- 
der the specified loading are obtained from considering the equilibrium of an element 
cut out by two diametrical sections and two cross sections perpendicular to the axis 
of the cylindrical shell as shown in Fig. 1. 

In this discussion the usual assumptions are made; namely, that the material is 
isotropic and follows Hooke’s law, the undeformed tube is cylindrical, the wall thick- 
ness is uniform and small compared to the radius, the deflections are small compared 
to this thickness so that second order strains can be neglected, and that straight lines 
in the cylinder wall and perpendicular to the middle surface remain straight after 


distortion. 
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The notation used for resultant forces and moments per unit length of wall section 
are indicated in Fig. 1. After simplification, the following equations of equilibrium are 


obtained :* 












































ON, ON¢:z OM.zg OM,y 
= 0, a — + as = 0, 
ax ag ax ag 
aN. aN x6 aM, Moz 
amunipias a — = 0, a _ g™= 0, (1) 
ab -— = ax a¢ a0 
00, 00 
a = + Rs + Not ga = 0, (N zg es Noz)a = 0, 
Ox Od 
in which g is the normal pressure on the element. 
If 0, and Q, are eliminated from Eqs. (1) and the relations 
Neo = Nosi Mag = — Mos 
are used, the six equations in (1) can be reduced to the following three: 
ON, ON 
Ox 0g j 
oN ON; OM, 1 dM 
*+a—*+—*-_— —*=0, (2) 
0g Ox Ox a 0d 
2 OMe aM, 1 0M, Ns 
— — —— 4 4 5 4g 
a 00x Ox? : a 


The relation between the resultant forces and moments and the strains of the 
middle surface will be taken the same as in the case of a flat plate: 





V Eh (e. + y Eh — y y vyEh 
Ns = i_ ya Ez VEw), Vg = 1 o y2 ey VErz), Veco = Voz = 21 rs y) 
M.= —D(X.+0X,), Me=—DXo+X.), Meoso= —Mee= Dl — Xen 


where D = Eh*/12(1—v?) is the flexural rigidity of the shell and h is the thickness. 

Resolving the displacement at an arbitrary point in the middle surface during de- 
formation into three components—ux along the generator, v along the tangent to the 
circular section, and w along the normal to the surface drawn inwards—one finds that 
the extensional strains and changes of curvature in the middle surface are 








Ou 1 dv w dv Ou 
a a cara ar 16 ax aah 
; 0*w ; 1 0’w 1 dv r 1 d°w 1 dv 
camry “o" ae 2 te en wares 


Hence, Eqs. (2) can be put int the form of three equations with three unknowns 


u,v, W: 





07u i+v 0% vy Ow 1—v 0% 
2 dsdx a Ox 2 os? 


* See Ref. 5, p. 440. 
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asx Zz Ox? a os 
(3) 


h? 0*w Ow h? 070 07n 
+—( — +—-) + (a-» +=) -0, 
12a \Ox?ds Os? 12a? Ox? Os? 


h? 1 /ov w Ou h? d%y 0%y 1 — vr? 
5 Vw -—(~-=4+ \+ (2-» + )-= —q=0, 


aaa a 
12 a \@s a Ox 12a Ox*ds Os? 




















where 
Ss = dad. 

In the problem under investigation, the quantities u and v are of the order of 
magnitude of \/hw/a. Consequently the last term and the third term in the second 
and the third equations, respectively, in (3) can be neglected safely. 

In order to solve the simultaneous equations (3), one can apply first the opera- 
tion 0?/dx?, and then 0?/ds? to the first Eq. (3). Solving in each case for the term con- 
taining v, and substituting these expressions in the equation obtained by applying 
0?/dxds to the second Eq. (3), one obtains an equation from which v has been elimi- 


nated: 











aV‘u4 = vp 


Ow Ow i+v kh?/ Ow Ow 
(fot) 


_—— — + 
Ox? Oxds? 1 — vp 12 \dx*ds? Oxds4 


Similarly, applying 0?/dx? and 02/ds? to the second Eq. (3) and solving for the terms 
containing u, and substituting in the first Eq. (3) after applying 0°/dxds to it, one 
obtains an equation from which u has been eliminated: 


(5) 














0*w aw ih? 2 Ow 3—v dw ew - 
aV‘y = (2+ ») : 


rs pat _ teas! 
0x7ds Os8 12\1 —v dx‘ds 1 — v Ox?ds* as 


Applying 0/dx to Eq. (4) and 0/ds to Eq. (5) and substituting these two equations 
into the third Eq. (3), after applying V‘ to it, one obtains an equation from which 


both uw and v are absent: 





4 2 =») tw Fits _ Os ise , oa ee oe 
839 - —$—<—_—— —_ + —[ - (2 + ») —— : v -— = 0. 
a*h? Ox! a* \ds® dxtas? 0x*ds* D . 
It is evident that the third term in Eq. (6) is neligible in comparison with the other 
terms. Equation (6) is reduced to 
12(1—»’) dtw 1 


ep ot (6a) 
a*h? ox* D . 





Equation (6a) differs from the differential equation of the flat plate only by the 
second term. The flat plate equation can be obtained from equation (6a) by the sub- 
stitution of a= «. Consequently, this second term represents the effect of curvature 
in the problem of the cylindrical shell. 

3. Infinitely long cylinder loaded with two equal and opposite forces. The above 
equation will now be applied to an infinitely long thin cylinder loaded, as shown in 
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Fig. 2, by two equal and opposite compressive forces P acting at the ends of a vertical 
diameter. 














et eae 
v — 


7) 


Fic. 2. Loads and components of displacements of an infinitely long cylinder. 


The difficulties of integrating Eq. (6a) for this type of loading can be circumvented 
by replacing the concentrated force P by a distributed load g expressed as function 
of the longitudinal and circumferential coordinates, and applied to a small area which 
subsequently is reduced to an infinitesimal. This is made possible by representing 
the function in the longitudinal direction, by a Fourier integral and in the circum- 
ferential direction by a Fourier series. Since g is an even function of both x and s, it 


can be expressed by 


g(x, s) = E + > Gn COS 21 f 109 nace. (7) 
2 a 0 a 


m=2,4++* 


The displacement w can be expanded in a similar manner in terms of a function 
w(A) as yet undetermined: 


3) 


i alia Ax 
w= > cos f w(d) cos — dx. (8) 
a 0 a 


n=0,2+°> 


It can be shown that the above expression for w satisfies the following requirements: 
at the point where the load is applied, the deflection and moment are continuous, and 
the slope of the deflection curve vanishes. Furthermore, the deflection vanishes at 
infinity. Substituting Eqs. (7) and (8) in the differential equation (6a) one obtains 
the following relations. For n=0, 


[fo [CY 3Q)]- Sm 0)} ans 
J 9 a aD\a 2D a oo 


(qo/2D) f(r) 
wr) => ° 
(A/a)4 + Eh/a?D 





therefore, 





Similarly for n=2,4,---, 


(qnf(d)/D) [(A/a)? + (n/a)?}? 





w(A) = —— 


Hence, the solution of Eq. (6a) is 
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1 ? gof(A) x 
w= — cos — dr 
2D Jo (A/a)* + (Eh/a*D) a 
rm 7s EAN) TG/OP ng 
D » 2 a n/a)? + (n/a)?|4 + (Eh/a2D)(d/a)* a 


It is next desired to find g, and f(A). In order to accomplish this the functions q, 
and f(A) must be determined from the loading condition. This is shown in Fig. 3. 
Since the cylinder is loaded symmetrically with respect to the generatrix and with re- 
spect to the circle passing through the origin, only the positive direction need be 











considered. 
ol 
=-Cco re) 
——Xx 
Fic. 3. 
From Eq. ( 
1 - x Ax x 
(= -) = r f(A) cos aM, f(A) = f q *) cos — (=), 
Tv — a a a 
and 
x x 
(=) =1 when —62524S6, (=) = (0 when x>déandx< —6 
a a 
Therefore 
2 6 


2f* dz x : 
fo) = = f cos d(—) = sina —. 
TJ 9 a a 7 a 


Similarly g, can be determined from the expansion of the loading function along the 
circumference in a Fourier series. With 
w/2 


2 
6.3 = q(z)dz, 


Tv —7/2 


where z=s/a, and if z=7/2, s=7a/2, one obtains 


2 “ 4c 2 : Ss 4q . c 
oe = qds = —4q, m=-—f g cos n — ds = — sin nm —» 
wa wa Tad _¢ a nN a 


where ¢ is as shown in Fig. 3. Substituting g, and f(A) in Eq. (9), one finds that 








1946] CYLINDRICAL SHELLS SUBJECTED TO CONCENTRATED LOADS 19 











1 *” (8gc/m7aX) sin \d/a Ax 
w= — —_—-——- —— cos — dd 
2D . 0 (A/a)4 + Eh/a*D a 
1 — ns (°* (8q/x*md) sin nc/a sin \5/a[(A/a)? + (n/a)?) Ax 
+- >> cos— cos — dx. 
D sihhes- adJy [(A/a)? + (n/a)?]4 + Eh/a*D(d/a)* a 


Next, the case of a concentrated load applied at the origin may be considered. 
Such a load can be obtained by making the lengths 26 and 2c of the loaded portion 
infinitely small. Substituting 

ee _ ne ne 


P = 4q¢6, sin — ~ —» sin — = — 
a a a a 


in the above equation, one obtains 





w= — 


v= 


Pa’? i” cos A(x/a)dd 
0 








DrJo M47? 
we ns (°° [d? + n?]? cos (Ax/a)dd 
cos — , (10) 
7D au2i--- G6 [\?2 + m?]4-+ Jr4 
where 
2 


Eha? a 
J? = —— = 12(1 — »’) (=). 
D h 


In order to evaluate the definite integrals in Eq. (10) Cauchy’s theorem of residues 
will be applied. Let us consider the integral 


f * cos A(x/a)dr 
————— ee 
0 rt + J* 
where the characteristic equation \4‘+J?=0 has four complex roots 


d= JU — 1) M4, 


Cauchy’s theorem yields 


= cos \(x/a)dr vy — ‘JT x E ‘J x 
f — = — = J-ire VI 12212) | cos — —-+sin — —/j. (11) 
J 9 4+ J? 2/2 26 2a 


The rational function in the integrand of second definite integral in Eq. (10) can 
be expressed in the form of partial fractions, 


(A? + n2)? 1 1 ‘ 1 
[Ae + m2]#+ sat 2 Ltt n~?+ ir 24 02)? - Wr? 











1 1 1 1 
 Qug Bag a? — a® 2? — a? 
1 2 1 2 1 
Fiteeaanb mae (12) 
a? — a? \? — a? buat Bad 
3 4 3 3 
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where +a;, tae, +a@3 and +a, are the roots of the denominator, 


t+a,= ta*#=+A+ iB 


rs ! 
=2l/cwtors(-F0)-o- 9] 


J \ 
lV iia +6)+o—], 
V2 2 
asx 





+ 


(13) 


i ) 


2 j 
ye soni +6) — (t+) 
‘ i 2 2 J 2 
| 4/ + 7) +(F+6)+O +0], 


where the asterisk denotes the complex conjugate, and 


|4 
R 
Il 
+ 





ll 
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|+ 























o= Vi(R?+ WV, n= Vi(R2 — J), R. = m2JV/1 + (J/4n?)*. (13a) 
Hence 
] fo (2 + n?)re2ledh 
JL (N+ 4 + TM 
2mi( ae 4) es tahatons 
= — {— -(n — ip eta /a oa —— (yn — ip e'taz/a 
st Qa1a2 a1 
a4 ate) a3 ee 
— —— (n+ ip)et*as/¢ + — (n + ip)e'*/¢ f- E (14) 
A3Q4 A304 ) 
Since 
ajag=- n? = 304, 


Eq. (14) can be simplified as follows. Now 





j © (2 + n?)? cos xd/adr 
0 (A? + n7)4 + J*r4 





Ax Ax 
= ize | coc + 7G) cos — + (¢G — 7C) sin =). Bra 


Cx _ Cz =e 
+ | (oa — B) cos — + (nA + $B) sin — | aay , (15) 
a a 


Simplifying the integrals (15) and (11) in Eq. (10), one obtains 


w/h VW3(1 — “(- ) ‘( /J x aul “f= “) File) 
——eemern cos 4/ — — + sin — —— fg VI Issle 
P/E 2x h eS a 2 a 


6(1 — v? 2 bed ) Ax 
Pe. (s ) > ent {[ wc + 16) cos— 


T nenG, 4-0 Ren? 
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, F Ax Cz 
+ (¢G — C) sin — | e-*/¢ + | (GA — 9B) cos — 
a a 


Cx]. 
+ (nA + $B) sin <*] ariel ; (16) 


It is seen that the first term of the above expression is very small as compared 
to the second term, and therefore can be neglected without appreciable error. For a 


40x 10? 





| MAXIMUM DEFLECTION PARAMETER _ 
35x10” OF CYLINDRICAL SHELL 
VS 


30x10” 
25x10? 
20x10? 


5 x10” 


g. 
h 


Fic. 4 


certain value of the a/h ratio, G is found to be very large as compared to B. The terms 
containing e~°*/* can then be completely neglected, provided that x/a is not near zero. 
In the case when x/a=0 Eq. (16) can be simplified as follows: 








w/h | 3/2(1 — (=) > cosns/a V1 +2 a7) 
P/Eh? 5 2;am0 7 T OF catinss | = 
where 
3(1 — v*)a? 
=? = 1 + ————__ 
4n*h? 


The left side of (17) has a maximum at s=0. Figure 4 shows the variation of this 
maximum with the ratio a/h. Figure 5 shows the variation of w along the generatrix 
through a point of loading, and Fig. 6 shows the projection of lines of constant w on 
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the plane through the axis of the cylinders and perpendicular to the line of action of 
the two forces P. 

4. A cylinder of finite length loaded with two equal and opposite forces. The ex- 
pression for the radial deflection in a thin cylinder of finite length can be obtained 
from Eq. (16) by using the method of images.* If one imagines the cylinder of finite 
length prolonged in both the positive and the negative x-directions, and loaded with 
a series of forces, P, of alternating sense, applied along the generatrix (s/a=0) ata 
distance / from one another (see Fig. 7), then the deflections of the infinite cylinder 


CURVES OF CIRCULAR CYLINDRICAL 
SHELLS ALONG THE GENERATRIX 
(INFINITE LENGTH) 
+0 


“4x0 
“2x0 


2x0 
4xi0° 
£ 
‘ef 6x0 
s|<£ 
8x0 
104% 
2x0 
4x0 
6x 
Bx 


Fic. 5 


are evidently equal to zero at a distance //2 from the applied loads P. Hence one may 
consider the given cylinder of length / and radius a as a portion of the infinitely long 
cylinder loaded as shown in Fig. 7. From Eq. (16) one finds that the deflection of 
any point, 6, (at a distance [ from the s-axis) on the shell due to the load P acting at 


the center is 








P. 2 © 
wo eS ERI CoC + 10) cos A= + (66 — 20) sin a = eave 
2rD m=2,4-++ Ron a a 


a | (oa — 7B) cos cS + (nA + $B) sin c=] easel : (18a) 
a a 


The deflection produced by two adjacent forces a distance / apart is 


* This method was used by A. N4dai, Z. angew. Math. Mech. 2, 1 (1922), and by M. T. Huber, 
Z. angew. Math. Mech. 6, 228 (1926). 
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Pa? 2 cos ns/a {[e +6) yp l-—¢ 
Ww,=> A 1) COS # 
° SoD okt... Be 


aut | 
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Sood j 
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l fe 
+ (¢G — nC) sin A = = e-aceouel ; (18b) 


a 








l 
+ (¢G — yC) sin A Jereon + | + 7G) cos A 





Since the terms containing e~°'+®/¢ are all small compared to the other terms, they 























Fic. 7. Series of equidistant opposite forces acting on an infinitely long cylinder. 


can be neglected without causing appreciable error. One obtains similarly w., wa, - + -. 
The total radial deflection at any point 6 is given by the sum 


W= Wat Wtwet::: 


Pa? = cosns/a c P § 
=— ) - {| 4 — nB) cosC — + (nA + $B) sinC =| e-Gtla 
2xD n=2,4++> Ron a a 





l 
+ (¢C + 7G) | cos A £ e-Btle — 2cosA Ls cosh B © (cos. — e Bila 


a a a a 


2Al l 
— cos —— e~*Bl/a 4. cos 3A — ec 3Blla — .... 
a a 
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a a a a 


; oss EP l R l 
— 2 sin A — sinh B = (sin A — e Bile — sin 2A — e*44/e 4... )| 


sae f S Pepa 
+ (¢G — 7C)| sin A — e-¥!/¢ — 2 cos A — cosh B —| sin A — e 84/4 


a a a a 


! l 
— sin 2A — e~*B!/e + sin 3A — e38l/a — ... ) 
d a 


oo l l 
+ 2 sin A — sinh B = (cos 4 —eBlia — cos 2A — oe *Blla 4... ) | : (19) 
a a a a 


We sum the series in the above expression, obtaining 


eo 


: * [e-m(B- iA)I/a 4 em (B+iA)l/a] 
m 


=1,3--- 


~ l 
> ele cos mA — = - 


m=1,3-+- a 


j= bdo] = 


sinh (Bl/a) cos (Al ‘a) 


—% x (Bla) cos? (Al/a) + cosh? (BI/a) sin? (Al/ xh 


l e~'Ble [sin (BI/a) cos? (Al/a) — cosh (BI/ a) sin? (4 Al/a) | 


> 6 ™842 cos mA — = - # 
ao a 2 (sin? (Al/a) cosh? (Bl/a) + cos? (Al/a) sinh? (BL ‘a)| 























2 l 
> eo ™Ble sin mA — = 
m=az1,3-°° a 


1 Cs) 
Pike > & [e-m(B- iA)l/a — ¢ m(B+iA)t/a] 
21 m=1,3--- 
1 


|: cosh (Bl/a) sin (Al/a) =| 
~ 2 Lsinh? (Bl/a) cos? (Al/a) + cosh? (BI/a) sin? (/ Al/ a) 
e844 cos (Al/a) sin 1 (4 Al/a) [sinh (Bl/ a) - + cosh (Bl, ‘a)] 


> ele sin mA — = eee 
wales a 2[sinh? (Bl/a) cos? (Al/a) + cosh? (BI/a) s sin? 2 (, Al/a)| 

















Thus Eq. (19) is reduced to 


w/h 6(1 — v?) ( a ) * COS NS/ /a j c 
= — — ( anit A 
P/Eh? 7 h -> .. Reon? if we + #6) 








+ (¢G — 7C) sin A =] gic 
a 


"y - 
+ | (os — 7B) me — + (nA + $B) sin c=] e-Gtla 


sinh (Bl/a) cos (Al/a) — e~8"/2[sinh (Bl/a) cos? (Al/a) — cosh (Bl ‘a) si sin? * (Al ‘a)| 





sinh? (Bl/a) cos? (Al/a) + cosh? (Bl/a) sin? (Al/a) 





4 | (6 — nC) sin A Lg sinh B a2 (@C + nG) cos A z. cosh B y 
a a a a 


cosh (Bl/a) sin (Al/a) — e~8"!* cos (Al/a) sin (Al/a) [sinh (Bl/a) + cosh (Bl/a) | 
sinh? (Bl/a) cos? (Al/a) + cosh? (Bl/a) sin? (Al/a) 





> 4 | (6 — nC) cos A £ cosh B £ + (¢C + 7G) sin A Ls sinh B =| . (20) 
a a a a 
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It is obvious that the first two terms of Eq. (20) are equivalent to the solution of 
the infinitely long cylinder given by Eq. (16). The remaining terms are evidently the 
correction factors due to the restrained edges at the two ends of the cylinder of finite 
length. The radial deflection under the applied force can be obtained by putting 


¢/a=0, 








w/h 6(1 — »*) fa\? & cosns/a 
= (+) A} (6C + 1G) + (A — 9B) 








(P/Eh?) in T ne2,s--- Ren? 
~ C + 16) sinh (Bl/a) cos (Al/a) — e~®”[sinh (Bi/a) cos® ‘Al/a) — cosh (Bl/a) sin? (Al/a)] 
; sinh? (Bl/a) cos? (Al/a) + cosh? (Bl/a) sin? (Al/a) 
“enw cosh (Bi/a) sin (Al/a) — e~Blle cos (Al/a) sin (Al/a) [sinh (Bl/a) + cosh (Bl/a)]} 
sinh? (Bl/a) cos? (Al/a) + cosh* (Bl/a) sin? (Al/a) 





sinh (Gl/a) cos (Cl/a) — e~@¥/ [sinh (Gl/a) cos? (Cl/a) — cosh (Gl/a) sin? (Cl/a)] 








— sinh? (Gl/a) cos* (Cl/a) + cosh? (Gl/a) sin® (Cl/a) 
eee cosh (Gl/a) sin (Cl/a) — ¢~Gle cos (Cl/a) sin (Cl/a) [sinh (Gl/a) + cosh (Di/a) } (21) 
stale sinh? (Gl/a) cos? (Cl/a) + cosh? (Gi/a) sin® (Cl/a) 


Some applications of the solution of the problem of the infinitely long cylinder. 
The problems of a couple acting on an infinitely long cylinder in the direction of 
either the generatrix or the circumference can be analyzed by using the solution given 
by Eq. (16) for a single load. The action of the couple is equivalent to that of the two 
forces P shown in Fig. 8, where limaz.»9 PAx = T,. 














——X % 


_—™, 
on f 4. 
+ 
+ 
| 





Fic. 8. Two couples acting on an infinitely long cylinder. 


It is easy to see that the deflection for the case when the force P is applied at the 
point QO, at a distance Ax from the origin, can be obtained from the deflection w, 
given in Eq. (16), by writing x —Ax instead of x and also —P instead of P. This and 
the original w are then added. The radial deflection due to the two equal and opposite 
forces applied at O and QO, is now obtained in the form 

— wr = w(x, s) — w(x — Ax, s). 
When Ax is very small, this approaches the value 
dw(x, s) 
——— Az. 


Wr = 
dx 


As T, is the moment of the applied torque and is equal to PAx, the radial deflection 


due to this torque is 
T. dw 


i sn cay (22) 
P dx 


WT, 


where w is the radial deflection due to the concentrated load P. 
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For the radial deflection due to the couple acting along the circumferential direc- 
tion one finds similarly (Fig. 8) that 


T. dw 
Fr a@s 


UT, = 


(23) 


Substituting w from Eq. (16) in Eqs. (22) and (23) one obtains for the couple act- 
ing along the circumferential direction, 





wr/h _ 6(1 — »*) a 2 2 cosns/a poma - : 
T/Eh He Rn? {e##/2 cos Ax/a[A($G — nC) — B($C + nG)] 


— ¢B2l0 sin (Ax/a)[(¢C + 1G)A + B(OG — 9C)] 
+ e-9l cos (Cx/a)[C(nA + $B) — G(GA — 9B)] 
— e2/@ sin (Cx/a)[C(¢A — nB) + G(nA + $B)]} (24) 


while for the couple acting along the generatrix direction, 


wr,/h 6(1—v%)/a\? & sinns/a 
= SEE) dR 6c + 16) 08 (42/0) 








T./Ek : Ve. oe 
+ (¢G — nC) sin (Ax/a) ]e~8*/* + [(6A — B) cos (Cx/a) 
+ (nA + @B) sin (Cx/a) ]e-@=/+}. (25) 
In the case when x/a=0, 
wr;/h 
= 0, atany s/a. 
T./Eh* 


Hence the condition that the slope of the deflection curve dw/dx must vanish under 
the concentrated load (x/a=0) is satisfied. 
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THERMAL DEFLECTIONS OF ANISOTROPIC THIN PLATES* 


BY 


WILLIAM H. PELL** 
Bell Aircraft Corporation 


1. Introduction. Equations governing the deflection of an isotropic thin elastic 
plate subjected to a temperature distribution of the form 


T(x, y, 2) = To(x, y) + 2T (x, y), (1.1) 


where the neutral plane of the plate is taken to lie in the xy-plane, have been derived 
by Nadai.' He did not consider the solution of these equations, and they have not 
been treated to any considerable extent by subsequent writers in thermo-elasticity.? 
The isothermal theory of anisotropic thin elastic plates has been developed principally 
by Boussinesq,* Voigt, and Lechnitzky.’ It appears that the only treatment of 
thermal effects for the anisotropic plate is due to Voigt,* who considers a simple case 
in which no bending of the plate occurs. 

The first part of this paper is concerned with the derivation of two partial differ- 
ential equations governing the deflection of a thin elastic plate possessing one plane 
of elastic symmetry parallel to the faces of the plate, and subjected to a temperature 
distribution described by a function of the form (1.1). One of these equations, with 
suitable boundary conditions, defines a stress function F; the other, the deflection 
function w. In the second part, recent results in anisotropic plate theory are used to 
solve the equation for the stress function with rather general boundary conditions 
for the case where To(x, y) isa polynomial in x and y. The problem of solving the equa- 
tion of the deflection is a difficult one, and a solution valid throughout the region 
enclosed by the plate is not available. Since the thermal deflection problem for the 
isotropic plate is of interest in itself, however, the case of the isotropic circular plate 
with radial temperature distribution is considered in the concluding portion, and the 
solution is obtained. 

2. The thermo-elastic equations for anisotropic plates. Let us consider a thin 
elastic plate composed of a medium possessing at each point at least one plane of 
elastic symmetry parallel to the middle plane of the plate, which is chosen to lie in 
the xy-plane. Let the plate be subjected in its interior to a temperature distribution 
given by (1.1). The plate is supposed acted upon by forces on its edge lying in the 
middle plane, but to be free of lateral load and body forces. 
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The Kirchhoff assumptions of the thin plate theory lead to the well-known funda- 
mental relations 








Ow Ow (2.1) 
“u=-—2 , = —3z- , 2. 
Ox oy 
| T22| K max ic! teak | rey | I, (2.2) 


valid throughout the thickness, 2h, of the plate, where w is the deflection of the middle 
surface of the plate, and u and v are the displacements of a point (x, y, z) of the plate 
in the x- and y-directions, respectively. To the assumptions (2.1) and (2.2) is adjoined 
the following one: the stress tensor 7 at any point in the plate is the sum 


r= 7°+ 7! (2.3) 


where 7° is a plane stress tensor generated by To(x, y) and reactions at the plate edge. 
and r' arises from the bending of the plate, i.e., from the action of z7;. It should be 
remarked that this supposition is fundamental in the thermo-elastic theory of thin 
plates presented here. Together with (2.1) and (2.2) it serves here the same purpose 
that (2.1) and (2.2) do alone in the isothermal theory, i.e., they reduce the thermo- 
elastic plate problem to one two-dimensional in character. 

The generalized Hooke’s law’ taking into account thermal effects is expressed by 





Ou 
— = Ay T22 + Qy2T yy + QisT2z + Q16T zy + aT, 
Ox 
Ov = 
he = Qi2T22 + Q22T yy + d23T22 + 267 zy + a2T, 
y 
Ow sig 
= Ay3T zz F Go3Tyy + 33722 + A367 zy + a3T, 
dz (2.4) 
Ov Ow : 
—— ot —— = Cutys + Q4sT 22 + 2a,T, 
Oz Oy 
Ow Ou re 
— aL — = W45Tyz + Q55T zz + 2asT, 
Ox Oz 
Ou ov = 
‘io + ae = Qi6T zz + Q26T yy + 1367 zz + 667 zy + 2aeT, 
oy Ox 


where the a;; are elastic constants and a; are the coefficients of thermal expansion of 
thermal expansion of the medium. For the plane stress system® (2.4) reduces to 


ou® 0 0 0 
ae = Q11Tzz + Qy2T yy + 2167 zy + aT», 
Ox 
dv® 0 0 0 7 
— = AyoTrz t+ G22Tyy + AeeT ry + a2T 4, (2. 5) 
oy 

ou ov® 0 0 0 

<2 ot = = AyeT ex oS 226T yy + (2667 zy + 2aeT 0. 

oy x 


7A. E. H. Love, Mathematical theory of elasticity, Cambridge University Press, Cambridge, ed. 4, 


1927, pp. 151-160. 
8 Displacements, strains, and stresses associated with 7» will be denoted by the superscript 0; those 


associated with 27;, by the superscript 1. 
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With the introduction of a stress function F(x, y), one has 


0 o°-F 0 0°F 


0 °F 
ton SF Tz ’ Tus = 
‘Oy? . Oxdy ax? 


The resultants of these stresses acting across the thickness of the plate are, respec- 


(2.6) 








tively, 


Neo 2htemn MNu@2hten Nea titer 


Since the displacements must satisfy a compatibility condition, it follows from (2.5) 


and (2.6) that F must satisfy the equation 


o*F ‘ o*F o*F 
26 antOy + (2412 + a6) ——— ax°dy 2 — 2ar6 axdy? 11 ays 
0°T 0°T» 0°T» 
AN \" ax? ‘a sf 











It is assumed that the derivatives appearing in (2.7) are continuous. 
The generalized Hooke’s law for the strains and stresses associated with the tem- 


perature function 27, permits one to write 


1 
zz 


T 


zy 


Ou} 
C11 —— a 12 
Ox 


1 Tou! 
+ ou(> [+ 
2 Loy 


dv! 
:) + éno( — aT, 


Ou} 
= €12| -—— — a3z7 
Ox 


4 f4 [ 4 
od ay 


ll 


du! 
tai — Cf 
Ox 


1 fou! 

+ C36 (> E + 
2Ldey 

7 dv! 

r:) + eao( — a2 

dy 
1 
| - eT), 


Ou} 
C164 7 — 232 
Ox 


1 [du 
Ot St > 
2 Loy 


The assumption (2.1) gives 


ou} 


Ox 


0*w dv! 


— & 


nw 





. pn 
Ox? oy 


Ov) 


Ox 


dv! 


Ox 


Ox 


Ov 


Ox 


; Ov 
T1) + ¢12| — — ao2T, 
dy 


| - aT), 
oy 


| - owl), 


/dv'! 
r:) a | a oe ao2T 


\ dy 


dot r ) 
ae | Oe 8 ee 
J 


0*w 1 


dw! 
) + Cis (= 
Oz 


Ow! 
) + C23 (= 
Oz 





dw! 
) + Css 
Oz 


dw! 
r:) + C36 (= 
dz 


ou} dv! 
ae See 
oy Ox 


_ assT1) 
_ axsT1) 
_ axsT1) 
_ axe) 


“axdy 





(2.7) 


(2.8) 


(2.9) 
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and (2.9) applied to the third equation of (2.8) yields 


Ow z 0*w 0? w 
i a3zT SS oe at Oi ee a,l oa C23 (= ed aol ) 
j UV" Nex? ay? 











02 C33 
0*w = 
+ C36 —_ agl Fig (2.10) 
Oxdy j 
Now (2.9) and (2.10). are inserted in (2.8) with the result 
1 0*w 0?w 2w B 
Ter = — 3( du or bis >} bis —— + aT }, 
Ox? Oy? Oxdy 
1 0*w 0*w 0? 
Ty = — S&S bie —- + bee nme sie Dag ee i gk 9 ’ (2 11) 
ox” Oy" Oxdy 
1 0*w O*w 0*w = 
Tag 7? = 2 bis — + bee ; + 566 gee -+ a6T }, 
Ox* ¢ Oxdy 
where 
CisCks ; 
bix ys =) 1, k _ t; 2, 6, 
C33 
aj = by:01 + doide + Ddeide, 1 = 1, 2, 6. 
Noting 7}, |.-+1=T7.|-+,=0, we may now integrate, with respect to z, the equa- 
tions of equilibrium 
na 1 1 f 4 = 
OT ze 4 OT sy OT as = 0, ts ‘ OT yy 4 OT ys aa 
Ox dy Oz Ox Oy Oz 
obtaining 
1 22 — h? Fw Ow Ow 
res = —| bu =, + Doe + (ie + bee) 
y 0x3 0x70 0x0" 
dw OT, OT, 
+ bx — + i + me ie 
Oy Ox Oy 
; (2.82) 
1 s*— }* d*w 0°41 Fw 
Tye = ———| bre — + (012 + Dee) — + Zhe —— 
2 0x8 Ox*dy Oxdy* 
0*w OT, OT; 
+ beg — + acs — + 22 — }- 
Oy? Ox oy 


The resultants of these stresses acting across the thickness of the plate are, respec- 


h h 
QO, -{ T2202, Oy -{ Ty. 
s =" 


h 


tively, 


The following equation expressing the condition of statical equilibrium of an arbi- 
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trary element of the plate may be obtained in the usual way’ since the derivation does 
not depend on the material composing the plate: 

00, dQ, 0*w 0*w 0*w 

=u ee “Ga” 2N zy —— + Ny — = 0. 

Ox oy Oxdy dy? 
The values of Q, and Q, obtained from (2.12) are now introduced, and the result 
is the following differential equation for the deflection: 








d4y 44 dw dw Ow 

bis = = 3b16 — “ts 2(b12 + ba) ——— + 3beg ——_ + bas — 
Ox 0x01 Ox*dy* dxdy4 dy* 
0°T, 0°T, 0°T, 

=— 1a) + 2as —* + a 
Ox? dxdy dy? 

W eile oe a (2.13) 
2h? ** ax? ‘  axdy Mia dy? 


Again, the continuity of the derivatives appearing is assumed. 

The portion of the xy-plane occupied by the middle plane of the plate will be 
called So; the boundary of So will be called Co, and it will be assumed that Cp is an 
analytic curve. The problem of thermo-elastic deflection is solved if a solution for 
each of (2.7) and (2.13) valid throughout S» can be found which satisfies appropriate 
boundary conditions on Cp. 

3. The stress function. Since derivatives of F appear in (2.13), the solution of 
(2.7) will be considered first. Lechnitzky’® has shown that the roots of the character- 


istic equation 
a;ys* — Zayqu® + (2a12 + Ge6)u? — 2Za2qu + dee = O (3.1) 


are necessarily complex. These roots will be denoted by nz=a% +78, R=1, 2, where 
6,0. Two cases must be distinguished: ype and pw, =pe. 

If F, denotes a particular solution, and mw; and pe: are distinct, then the most gen- 
eral solution of (2.7) is given by 


F(x, y) = Fy(z1) + Fi(é1) + Fo(z2) + Fe(Z2) + Fp(x, ), (3.2a) 


where F;(z,;) is in each case an arbitrary function of z,=x+pxy. The F;(z,) is analytic 
in the region S; of the z,-plane which corresponds to the region So of the z-plane under 


the transformation™ 


= piso t Gil, k = 1,2, (3.3) 


where 


pe=H1—im), g=HL—im), k= 1,2. (3.4) 


4 


If ui=pe, then 


® A. NAdai, loc. cit., pp. 233-235. 


10S. G. Lechnitzky, loc. cit. 
1 A summary of anisotropic plate theory in English can be found in I. S. Sokolnikoff’s Mathematical 


theory of elasticity, mimeograph lecture notes, Brown University, Providence, R. I., 1941, pp. 319-329. 
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F(x, y) = 2,F (21) + 21F 1(21) + G,(21) + Gi (21) +- F(x, y). (3.2b) 


For the plate without holes, S; is a simply-connected domain, and F;(z,) and G,(2;) 
are single-valued analytic functions in S;, and Sj, respectively. 
4. The boundary conditions for F. The plane stresses must satisfy the well- 


known conditions 
0 0 
= Tz, COS (x, m) + Tz, Cos (nm, ¥), 


X 
Vy 


(4.1) 


0 
n 
0 0 0 

n = Tzy COS (x, m) + Ty, cos (y, n), 


where Xz, Y, are the x- and y-components, respectively, of the force acting on the 
edge of the plate, and m is the exterior unit normal to Cp. If F is introduced through 
(2.6), the boundary conditions on F may be written 

oF _ OF 


. — 
Ox Oy 


= if (X, + i¥,)ds + ¢ = fils) + ife(s) + ¢, (4.2) 


where s is the arc-length along Co, measured from an arbitrary point with the usual 
convention as to positive s, and c=c’+ic’’ is an arbitrary, complex constant. A 
familiar alternative form of (4.2) is 

OF 


f(s), | (4.3a) 


on on Co, 


F = g(s), ) (4. 3b) 


where f and g are prescribed functions along Co, except for the arbitrary constants 
c’ and c’’ appearing in them. 

The solutions considered here are assumed to be such that Ff (z;), Fi(2:) and 
G{ (2) are continuous in So+Co. In this case (4.2) and (4.3) hold. To ensure the exist- 
ence of such a solution, it is necessary to demand that 


f X.ds = 0, f Yds = 0, (4. 4a) 
Co Cc 


{ fi(s) cos (x, s) + fo(s) cos (y, s)}ds = 0. (4. 4b) 
Co 


The physical significance of these conditions is of interest. Equation (4.4a) expresses 
mathematically the fact that the resultant of the external forces acting on the plate 
must vanish, and (4.4b) that the resultant of the external moments must vanish. 

The analytical similarity between the boundary value problem presented by (2.7) 
and (4.3) and that of the clamped plate under lateral load makes it possible to use 
recent results on the lattter problem obtained by Morkovin.'? His treatment depends 
essentially on the handling of the boundary conditions. 

The solution (3.2a) leads to boundary conditions expressed in terms of both 2, 
and 2, along either C; or C2 (the boundaries of S; and Se, respectively). Boundary 
conditions in terms of a single variable are obtained by mapping conformally a band 
of the 2,-plane containing C; in its interior onto an annular region of a {,-plane con- 


2 VY. Morkovin, Quart. Appl. Math. 1, 116-129 (1943). 
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taining in its interior the circumference of the unit circle y, in a way such that C; 
is mapped into y;. By a proper choice™ of the functions effecting this mapping, the 
transforms zg; and 2 of any given value of z on Cy can be mapped onto 7; and 72 in 
such a way that {:={2. This common value is denoted by ¢ =e**. Thus the boundary 
conditions contain ¢ alone. 


Let 
Se =or(Se), k= 1,2, (4.5) 
be the functions achieving the desired mapping. Then F has one of the forms 
F = o1(§1) + G1(F1) + 2(S2) + b2(F2) + Fo, (4. 6a) 
F = oy(Fi)bilSi) + @1(1) bE 1) + PilSs) + VilF) + Fo, (4.6b) 


where :({%) =F(we(x%)) and Wilf1) =Gil(wi(f1)). These functions are analytic and 
single-valued in some neighborhoods of yx, and hence possess Laurent expansions 


00 


ox(Sx) = Do Ynrtk Wilf) = > pasha (4.7) 


t-_—o n=—oo 
The coefficients y,, and uw»; are to be determined from the boundary conditions 


on F, expressed now on ¥; or ¥2, as desired. If one defines 


H,(c) Pwr (o)o + Gicwx (a)a, 


@e(o) (4.8) 
(pepe + ede) + ¢——— ol (o) 2piges 
then boundary conditions on 7; or Y2 (k=1 or 2, respectively) equivalent to (4.3) are 
given by 


d s OF 1 OF OF a OF oF 
of 8) faa (2- 12) ema(Z 42) | aon 
d 2B, Ox Oy / cy oy 


do on 


aF aF aF » , 
= a= sf | alo) (— —i >) — H,(@) (— +i =). ; +c’, (4.9b) 
oy Co Oy o 


where™ the left-hand member of (4.9a) is defined by 


_ ds OF 1 OF ae 
ww, —|.- — {yler) = + J (Gx) a (4. 10) 
doz, Onc, Bx don Co IoK Icy 


and c’”’ is an arbitrary real constant. 
5. The determination of F. The function F will be determined under the assump- 


tion that the edge force vector (4.2) has the form™ 


OF OF ~ 
(— + i) ee (5.1) 
2 7Co 


j=—ky 


J x(¢) 


1% V. Morkovin, loc. cit. 
\“ After the indicated differentiation of w,(o,), ¢,in F and its derivative are to be replaced by the com- 


mon value ¢ of o; and a2. Note also that ¢ =o. 
& Since ¢ =e", it follows that (@F/dx)c, and (@F/dy)c, are trigonometric polynomials of order ky. 
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where the c; are complex constants. The temperature function 79(x, y) is taken to be 
the polynomial 
k i 
0 j i-j 
Tis y= Ddt-wy (5.2) 
im0 jm0 
where the #;,;_; are real constants, and & is an arbitrary integer greater than or equal 
to 0. Then a particular integral 
kt+2 @ 
Fy = OD hi s—sriyi (5.3) 
t=2 jm 
of (2.7) can be found without difficulty. The Laurent expansions of w;(¢;), along with 
(5.1) and the ¢,-plane transform" of F,]c,, are now substituted in (4.9). The resulting 
equations yield recurrence formulas for the coefficients Ynx% (or Yn1 and fq1). 
For the case of the circular plate, the contour C> is a circle and the mapping func- 
tions are easily found to be 


tig 
Ze = wer) = on (t+), Yy,=—) k = 1, 2. (5.4) 
k 


Moreover, if wx({%) are given by (5.4), then a consequence of the single-valuedness of 
F,(z,) and G,(g:) is that (4.7) may be written 


r) 5 T, ee) " r. 
bx (fx) = Dru(si+ =), (01) = & saa +=). (5.5) 
nen : n=0 se 
The direct mapping function corresponding to (5.4) is 
Z2= ao (S.6) 
and hence it follows that 
k+2 a 
Frle, = dD. (Cro* + Cnr") = C(o). (5.7) 
n=0 


For k=0 the operator (4.10) becomes o(0/d0)+4(0/dc¢), and one then obtains 


ds OF , et aie 
io — | = > (D,o" + D,o") = D(c). (5.8) 
do OnJc, nm 
The coefficients C, and D, are in general complex constants. 
If A(c) and B(c) denote the right-hand members of (4.9a) and (4.9b) respectively, 
then the insertion of (5.1) and (5.6) in A and B yields 
qa *! ih a ‘1 a 
A(c) = - Dd (Ano® + Ane"), B(s) = = > (Bro* + Bie"), (5.9) 
n=0 n=0 


where 


16 This transform is most easily found by using the mapping function z=wo(fo), obtained by setting 
k=0 in (4.5), and defining uo=i. In this way one achieves a direct mapping from the z- to the {o-plane, 
and maps any point of Cy into a point oo with the same coordinates as that given by the successive map- 
pings (3.3) and (4.5). 
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( cr, n = 0, 
himtinhiains oF Se 
| dese O Rinks n= 2,3,---,h— 1; (5.10) 
> Cyt Cj | fel = «0, 
jak, J — ” 
B,=4 Go — até, sales 
Conti — Cat1 


—y n=2,3,---,k—1 





nN 


and the prime on >> indicates that the term corresponding to j =1 is to be omitted. 
For convenience in writing, let 


ke ko 
II(c) = >> (no + Tne"), A(o) = D> (Ano* + Ano"), ke = max [k+ 2, &, — 1], 


n=0 n=0 


where 
I, = A, — D,, A, = B, —C,, (5.11) 
with A,, B,=0 if n2k:, and C,, D,=0 if n2=k+3. 
If X28 and Y? are such that (5.1) is valid, then condition (4.4a) is satisfied, and 
(4.4b) demands that ¢ =. 
The determination of ¢; and y is simplified by using not the boundary conditions 
(4.9), but an equivalent set.!’? For the case of equal roots these are 


@(5)¢(o) + w(c) (4) + ¥(c) + Y(e) = ACO), (S.12a) 
w(a)o¢'(c) + w'(c)od(a) + o’(c) 


te 2 a 
“—— aS _ Ha) } {T@)a¥(6) ~ 3 cA'(c)a ~@I@} : (5.12b) 





Subscripts have been omitted in the above, since the distinction between the 2:- and 
Z:-planes is not involved in the discussion. Substitution in (5.12) of the series’® (6.4) 
for ¢(c) and ¥(c), together with w(c) from (5.4), then yields recurrence relations for 
Yn and up. For the explicit form of these rather lengthy formulas, the reader is re- 
ferred to the author’s thesis.!® 

It should be noted that equations (5.12) are valid for any given boundary of the 
admissible class, provided the functions II(¢) and A(c) corresponding to this boundary 
are found. 

The treatment of the case “pe does not depart materially from that of the case 
of equal roots, and hence will be omitted here.”° 

6. Extent of arbitrariness in solutions. Since a certain amount of arbitrariness 
is present in the functions occurring in the MuscheliSvili solution of the plane stress 
problem for isotropic media, it is reasonable to ask if this phenomenon persists in the 


17 V, Morkovin, loc. cit. 

18 That ¢and y may be so written is clear after reading Section 6. 

19 Thermal deflection of anisotropic thin plates, University of Wisconsin, 1943. 
20 For details see the reference of footnote 19, 
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case of anisotropic media. This question is answered affirmatively below, and the ex- 
tent of the arbitrariness determined for the case of equal roots. This arbitrariness 
stems in part from the fact that the stress function is here, as in the Muscheli$vili 
theory, the real part of an analytic function. In addition, it will be noted that the 
boundary conditions (5.12) contain three arbitrary constants, and one may expect 
this arbitrariness to be manifested in the solution. 

If the prescribed functions X%(s) and Ye(s) satisfy the conditions (4.4) and are 
representable in the form (5.1), then the uniqueness theorem”! for the first boundary 
value problem of elasticity assures one that for a given distribution of temperature T9, 
the state of stress in the interior of the plate is determinate. Since the stresses are 
given by (2.6), it is thus seen that the second derivatives of F are determined in Sp. 

If one lets F; be a solution of (2.7) and (4.3), lets Ff” =2Re { 2:F{” (2:1) +G{"(z:) } be 
another solution having the same second derivatives as Fy, and lets 


(2) ) (2 (2 
Fy’ =Fy—F,. = 2Re {aiF; (1) +Gi (2)} (6.1) 


then it follows that 
2_ (2) 2_ (2 2_ (2) 
OF, oF, OF; 
ax? Oxdy dy? 
throughout So. The real and imaginary parts of F? and G? must satisfy the Cauchy- 
Riemann equations in S,, and these together with the three equations (6.2) enable 


one to show that 


=0 (6.2) 





(2 - P ° 
Fy (21) = — insti + (m+ ine), 
Gy (21) (vs + ivg)zi + (v1 + ive), 


where v; and 7; are arbitrary real constants. Thus the functions (4.7) and 


r 
$1($1) = m + ine — insaps( ts T =) + 1 (f1), 
: (6.3) 


¥1(61) 


describe the same state of stress in So. The arbitrariness in (6.3) is removed by choos- 
ing v; and 7; so as to simplify ¢; and y:. The choice made here is such that 


Tr; (1) 
1 


Vi + ive + (vs + indapa (ts + F )+ Vi (f1), 


- © fC ™ t. 7 r 
$i(f1) = D(x + =) vit) = D pe (x + =) (6.4) 
n=] 1 n=2 , 


where Yu =Yu.- 
For the case 412, a more lengthy consideration®* shows that one may write 
n 


co — 2 ee 
oi(f1) = Dmx + =) $2(f2) = Dt (# + —), (6.5) 


n n 
n=2 1 n=l - 


where 22 = Y22. 
In both cases the arbitrariness is a reflection of that inherent in the values f(s) 





21 For details see the reference of footnote 19, 
% See the reference of footnote 19, 
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and g(s) which F must assume along Co. The choice of v; and ; implied by (6.4) (or 
(6.5)) is found to dictate the selection of the arbitrary constants in (4.3), or in (4.9), 
the modified form of these boundary conditions, and conversely. This somewhat in- 
verted method of eliminating the arbitrariness in ¢; and y; (or $1 and @¢2) is adopted in 
order to have these functions assume the form usual in the Muscheli8vili theory.?* 
7. The differential equation for the deflection; the associated boundary condi- 
tions. The coefficients N.:, Nyy, and Nz, in (2.13) may now be regarded as known, 
and the thermo-elastic deflection problem for the anisotropic thin plate is reduced to 
that of solving the partial differential equation (2.13) for w, subject to the appropriate 
boundary conditions to be satisfied on Co. For the first and second boundary value 
problems of plate theory these are, respectively, 
oH 





M, = m(s), 0, + P “= p(s), on Co, (7. 1a) 
s 
and 
Ow 
w= w(s), —=~w,(s), on Co, (7.1b) 
on 


where m, p, w, and w, are prescribed functions along Co. The quantities** M,, H,,. and 
Q,, are the flexural couple, torsional couple, and the shearing force, respectively, which 
act on the edge of the plate. 

The specialization of (2.7) and (2.13) to the isotropic case will now be given, since 
the resulting equations will be used in the sequel. For isotropic media 


1 2(1 + o) 
aa =—, ¢ I, 2, 3, ays = ——, 1=4,5,6, 
E E me 
ce «4. ot ps a (7.2) 
=—-—, 14j=1,2,3,4# J, a;;=0,1,7 = 4,5,6,1 #7, 


a; = a, t = I, 2, 3, a; = 0, ¢ = 4, 5, 6, 


where E denotes Young’s modulus, ¢ is Poisson’s ratio, and @ is the coefficient of 
thermal expansion. One easily finds that 











E E 
bu = bx = — —- bes = , 
1 — o? l+o 
Eo - 
bie TE Hee 9 bis => bag —- 0, (7.3) 
1 —o° 
ak 
Qa; > ara=— —9 ag = 0. 
1-—o 


Using these values of the elastic parameters, one obtains from (2.7) and (2.13) the 


equations 








VF = — aEy’*T», (7.4) 
3 (+ ov + — {a ae —\ (7.5) 
Vwur=- v*T — Nez 2N; Nw =: Se 
a o 1 D x2 y andy vy ay? 


23 1. S. Sokolnikoff, loc. cit., pp. 243-251. 
24 These are linear functions of M,, M,and Hzy. See I. S. Sokolnikoff, loc. cit., p. 326. 
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where D = 2Eh*/3(1—o?). Except for changes in notation, these are the same equa- 
tions that Nadai® obtains. The quantities (7.2) and (7.3) must also be inserted in 
(7.1) in order to obtain the boundary conditions for the isotropic plate. 

A solution of the deflection equation (2.13) valid throughout the domain So and 
satisfying the given boundary conditions on Co, is not available; the same is true of 
(7.5) and its associated boundary conditions.** By further specialization, however, it 
is possible to obtain the solution of the thermo-elastic problem for a case which is of 
some interest. This will be done in the following sections. 

8. The istropic circular plate with radial temperature distribution. Let us con- 
sider an isotropic circular plate of radius a, and let it be subjected to a temperature 
distribution given by 
T(x, y, 2) = To(r) + 27,(r), (8.1) 
where r= x+y’, and the origin is assumed to be at the center of the plate. As in 
section 2, the second derivatives of JT) and 7; are assumed continuous on 0Sr7 Sa. 
The edge of the plate is taken to be subjected to a uniform force P per unit length of 
the arc parameter s. 

It has been seen that before the deflection w can be found, it is necessary to solve 
(7.4) for F. In the case at hand, F, can be found easily whether or not radial symmetry 
exists, for if F is a solution of 

v¥ = — aET, (8.2) 
then it is also a solution of (7.4). But (8.2) is the well-known Poisson’s equation, and a 
solution is at once available from potential theory. Since T» has radial symmetry in 
the present case, however, the Laplacian operator becomes 


1d/d 
v= — “(r=), (8.3) 


dr 


and successive integrations of (8.2) give 
"dé cé 
F,(r) = - ak | =f T o(x)xdx (8.4) 
0 ¢ 0 


as the particular solution of (7.4) which is needed in F (see (3.2)). For the isotropic 
plate, the roots of the characteristic equation (3.1) are 4; =2 =7 and their conjugates. 
The transformation (4.5) becomes 


z= at (8.5) 
and if |¢| =p, then r=ap. If one lets f,(p) =F,(ap), then from (5.7) 
C(c) = f,(Voe) = f,(1), (8.6) 
and 
ds OF re) 0 
D(a) = to — —] = (« a as a=) 140) | = Voe fy (Voc) = fy (1), (8.7) 
Oa On Jc, Oo 0a ae 


% A. N4dai, loc. cit., pp. 264-268. 
26 Recent work by S. Bergman gives promise of extension to equations of the type (2.13) and (7.5). 
Bergman considers an equation a special case of which is 
Viu + atizr + 2bucy + C¥yy + duz + eu, + fu = 0, 
where a, 6,--- , f are analytic functions of x and y. See Duke Math. J. 11, 617-649 (1944). 
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i.e., C(o) and D(c) reduce to constants on ¥. This result is clearly a consequence of the 
fact that J» is a function of r alone. One finds that 
Co = 3f,(1), C; = 0, j=1,2,--: ’ 
Do = 3f7 (1), D; = 0, j=1,2,---. 


With the edge forces as specified above, it follows that one may write 


(8.8) 


0 = +6 
Xn+ iY, = — Pe 
and then (4.2) yields 
OF _ OF 
— +i1—) =aP(1—o)+ 6. 
Ox Oy / cy 
This is of the assumed form (5.1), and the constants defined in that expression have 


the values 


co = aP, co, = — aP. (8.9) 


These are now inserted in (5.10) and the result together with (8.8) substituted in 
(5.11). Proceeding as indicated in section 5, one finds that the recurrence formulas 


simplify to 
a 1 
n= -<{p 4 fy, yn = 0,n =2,3,---, un = 0, nm = 2, 3,---. (8.10) 
a 
The above results in conjunction with (3.2b) enable one to write 
1 cE ff’ : df i 
F(r) = ={- P+- -f To(x)xdx pr? — ak [ =f To(x)xdx. (8.11) 
2 r? Jo 0 Evo 


In view of the radial symmetry, it is expedient to write (7.5) in terms of polar coordi- 


nates in the form 


or\r 00 
1 dw 1 0°w 
+ No(— +S), oe. 


r Or r 302 


1 0?*w _ 0/1 dw 
iw = —a(l + o)V°7T, + = {¥en or + 2Nw — (- ) 


0 0 
where N,, =2hr°,, Nie =2hr*,, Nop =2h7?,, and 


0 1 a 1 oF ‘ es = . OF tii 
tre = — — Ht r= —-—l(|— — T= —- : 
, or oe 0 om ar\r 20 mo 


With F of the form (8.11), these stresses become 


0 ak a ak r 
Tr = — P+ —f To(x)xdx — ~—f To(x)xdx, 
a? r? Jo 


0 


I] 


Il 


ak f* ait f* 
eat die —f To(x)xdx + —f To(x)xdx — aE£T (1), (8.14) 
a 0 0 


0. 


in) 
4 
Il 
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Since 7}, N,,-, N+», and No are functions of r alone, V? has the form (8.3), and the equa- 
tion for the deflection may therefore be written 


d df1<d( dw d dT, 
ar drLr dr dr dr dr 
1 “[(? ak fre \xd \ 

— r—-— , uetied 

D ‘ r 0 oe dr? 


ak r dw 
+ (?. “Pe = =f T o(x)xdx —akE rs) es | (8. 15) 
- 0 


ar 


where P» is the constant 
ak 
Py) = — ” Tolx)xdx — P. (8.16) 
a? 0 
Equation (8.15) may be integrated with respect to r and an equation obtained 
thereby which is not only of lower order than (8.15), but which also has simpler co- 


efficients. Noting that 


WA? wae =f Tox)xdx) = 


ak : d*w ak _, \ aw 
= (Pe _- — Ta(x)xdx) — + (?. “4. —f T o(x)xdx — aET .) —) 
r dr? r? Jo 1 


0 ar 


we can carry out the desired integration immediately, obtaining 


kd Sen (- zl 
dr\r dr dr 
Oss 2h ak r dw ky 
= =— wl + 6) — “ (7. is a Tu(x)xdx) pone yr aoe (8.17) 
dr ie r? Jo dr r 
where k; is a constant of integration. Thus for the thin circular plate under uniform 
compression on its edge, and with T(r) and 7;(r) arbitrary save for certain conditions 
of continuity, the problem of finding the thermo-elastic deflection is reduced to that 
of solving the third order differential equation (8.17) with appropriate boundary con- 
ditions. 
A repetition of the above integration is impossible for T) and 7, of the general 
nature assumed above. Therefore, in order to complete the integration of the equation, 
T, and 7; will be supposed to have the form?’ 





m 1 n 
T(r) = to jr? Ti(r) = — ———_ DL’, (8.18) 
o(r) = =? ' a(1-+ oc) 2X " 


where ¢o;, 4;; are real constants and m, n are arbitrary positive integers. The solution 
of (8.17) may now be sought in the form of a power series. The polynomials 7 and 
T; are now inserted in the differential equation, and if one makes the abbreviations 


27 These polynomials may be regarded as approximations to power series representations of Ty) and 7. 
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bo; : 
harry dj = jhij, 





the result is 


dw 1 d*w m 2hP 1) dw » 
s+ + {Douay t = Dag, (8.19) 


dr® ry dr D r?) dr j=l 





j=0 


The obvious replacement 


dw 
u(r) = — (8.20) 


gives an equation of order two, and the further substitution 


2hP om 
bf =—  — Di dai 
D jl 


then yields 


1*4 1 du _ 1 a 
ees OF +(Lowi+ os - >) a= yo dsyari (8.21) 
j=0 


dr? r dr jul i 


as the equation to be solved for u. It will be observed that this equation has a regular 
singular point at r=0, and that the indices relative to this singular point are +1. 
Only the solution of the form 
oe 
u=r), dri (8.22) 
t=0 
will be considered here. This represents the solution relative to r=0 which is bounded 
there. It is evident from physical considerations that this boundedness must obtain 
for the simply-connected plate, and hence it is sufficient to consider the solution of 
the above form. The series (8.22) is now inserted in (8.21) and the following recurrence 


relation defining the A; is obtained: 


= due, J = —1,0,---,2 — 2, 
(GHAG H+ rset IA; + Do Aji; = { oe (8.23) 
i=0 0, jon—i,n,:-:- 


where A; =0, if 7<0. These equations permit \; to be expressed in terms of the arbi- 
trary quantity \» and the known quantities b¢ , b;, and d;. It is not expedient to give a 
formula for \;, since such an expression would be quite lengthy. The first several \; 
may be computed with little trouble if m and m are not large, but the labor involved 
increases rapidly as these numbers become larger. 

Up to now, the series defined by (8.23) is a formal solution of the differential 
equation (8.21). It is not difficult to show that this series converges for 0 Sr Sa and 
hence is actually a solution of the equation. For the discussion of the convergnece 
it suffices to consider the series 

oo 
> Aw, (8.24) 


tah 


where h=m-+n-2. If one defines 
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M = max {1, | bd |, | dil, | bel,---, | 5, | i, 
K = max { | dol, [ail,---, | Ana | I, 
and sets E(@) =(h+6)(h+2+86), repeated use of the recurrence relation (8.23) with 
j>n-—z2 enables one to deduce that 
K{M(m+1)} [Monee] 


: = Hk (m+2)+i 8.25 
E(0) E(m+2) E(2m+4) - - - E(k(m+2)+i) Mk (m+2)+ ( ) 





| Angicm+2)44 | 


for k=0,1, 2, - - - , 4 assuming the values 0, 1, - - - , m+1 for each k. It may then be 
shown that the dominant series defined by the picm42)44 converges uniformly for any 
finite r, and application of the well-known Weierstrass theorem yields the same con- 
vergence of (8.24) and (8.22). The desired solution w of the deflection equation is then 
obtained by inserting (8.22) in (8.20) and integrating. We obtain 


w(r) = 2? > xrit« (8. 26) 


t=0 
where x is a constant of integration and x;=\;/(i+2). The uniform convergence of 
this integrated series on the interval 0 <r Sa follows immediately from that of (8.22). 
Reference to the recursion formula (8.23) reveals that one may write 


ki = Aoki + 4; (8.27) 


where £; contains the parameters P, D, h, a, and some or all of the to;, while 6; con- 
tains not only these but also some or all of the #,;. Thus the deflection may be written 


in the form 


w(r) = Aowo(r) + wi(r) + «x, (8.28) 
where 
wor) = Déri, wilt) = Dbz’. (8.29) 
j=0 j=0 


It is clear that w;(r) is a particular solution of (8.21) and wo(r) the solution of the 
homogeneous equation which is bounded at r=0. 

9. The circular plate with clamped and simply supported edge. The constants 
Xo and « in (8.28) will be determined by the mode of support of the plate. The assump- 
tion of radially symmetric deflection limits one to the consideration of boundary con- 
ditions compatible with this type of deflection. The two methods of support most 
commonly encountered are those of the clamped edge and the simply supported edge. 
The boundary conditions for these are, respectively, 





dw 
=] = 0, W | rma = 0, (9. 1) 
a7 Irma 
and 
d*w o dw 
M,| rma = | + iat =| = 0, W | rma = 0. (9.2) 
dr? 7 7 |e 


The substitution of (8.28) in (9.1) and (9.2) yields for the plate with clamped edge 
the constants 
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awi’ (a) + ows (a) _ awit’ (a) + ows (a) 


awi"(a) + ews (a)’ “~ awi"(a)+ ews (a) MO 





A= - 








and for the plate with a simply supported edge, 


wy (a) wi (a) 
Wo (a) ihe we (a) Wo(a) —_ w,(a). 





These quantities are now inserted in (8.28), giving for the plate with clamped edge 
the deflection 
wy (a) 


Wo (a 


{wo(a) — wo(r)} + wi(r) — wi(a), (9. 3) 





wr) = 


and for the plate with simply supported edge 


thy ee tate ~ edt + ees hk. (9.4) 


aw)’ (a) + ow? (a) ' 

It will be assumed that w¢ (a) and awg’ (a)+ow¢ (a) do not vanish. It may be 
shown that an a which causes the former (latter) to vanish is the radius of the 
clamped (simply supported) plate with given temperature distribution T(r) for which 
P isa critical (i.e., buckling) load. The question of stability is not under consideration 
here, hence the above assumption is made. 

The case m=0 is of particular interest, for then w contains the Bessel function 
of the first kind of order zero. In this case T(r) reduces to a constant, and (8.21) be- 


comes 
d*u 1 du 2hP 1 n—1 
—-+— + ( —-—— — =) “u= } si d jar’. 


dr? Tr dr D r? j=0 





If the recurrence relation is written in terms of the «x; rather than the \;, the result is 


I 


: 2hP hue *#= —1,0,---,n— 2, 
(¢ + 4)*i+2 + ——K = ‘ . 
0 1 


D n—I1,n,---, 


from which one obtains easily the deflection 


2hP \ 
w(r) = Ao (4/— + wio(r) + x. 
BD Fs 


Here wyo(r) designates the form which w,(r) assumes for m =0. Proceeding as in the 
more general case, for the plate with clamped edge we find that the deflection is 


— 


w30(a) f [2hP 2hP 
w(r) = Jo ny ——) — Jol a4/ — 
2hP 2nP\ \ D D 
ae Jo a, a 
D D 
+ W30(7) aa W10(@), 


and for the plate with simply supported edge, 
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w(r) = — _awi' (a) + ows (a) 


J2HP ( /2hP 2hP /2hP 
hy > 20(24/ ) = a= (a4/—)} 
D \ D D D 
SJ =) ; ( /—) eee 
13 r4/ D Jo\a D + Wi0(7) W10(a). 


10. Conclusion. The problem of flexure is considered for a thin anisotropic elastic 
plate, subjected in its interior to a temperature distribution of the form 


T(x, y, 2) = To(x, y) + 27i(x, 9). 


The usual assumptions of thin plate theory, together with Hooke’s law extended to 
encompass thermal effects, permit one to derive two partial differential equations 
governing the deflection of the plate. For the isotropic plate these equations specialize 
to those given by Nadai. If thermal effects are supposed absent, and N.z, Nyy, Nzy are 
interpreted as arising from edge forces alone, then (2.13) becomes the equation for 
the deflection of an anisotropic thin plate stressed in its own plane. 

A method is given for determining the stress function F for a plate with edge 
forces representable in the form of a trigonometric polynomial, and the determination 
of F is carried out for the circular plate. 

For the thermo-elastic problem formulated with the above generality, a suitable 
solution of the deflection equation is not available. Accordingly, the problem is spe- 
cialized to the simpler case of the isotropic circular plate with radially symmetric 
temperature distribution and T(r), T7:(r) in the form of polynomials. The solution of 
the resulting deflection equation may be found in the form of a power series, and con- 
vergence established. Boundary conditions for the clamped and simply supported 
plate are then considered, and the deflection for each of these modes of support de- 








termined. 
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CONTRIBUTIONS TO THE PROBLEM OF APPROXIMATION 
OF EQUIDISTANT DATA BY ANALYTIC FUNCTIONS* 


PART A.—ON THE PROBLEM OF SMOOTHING OR GRADUATION. 
A FIRST CLASS OF ANALYTIC APPROXIMATION FORMULAE 


BY 
I. J. SCHOENBERG 


University of Pennsylvania and Ballistic Research Laboratories, Aberdeen Proving Ground 


Introduction. Let there be given a sequence of ordinates 
{yn} (n=0,+1,+2,-:-), 


corresponding to all integral values of the variable x =n. If these ordinates are the 
values of a known analytic function F(x), then the problem of interpolation between 
these ordinates has an obvious and precise meaning: we are required to compute 
intermediate values F(x) to the same accuracy to which the ordinates are known. 
Undoubtedly, the most convenient tool for the solution of this problem is the poly- 
nomial central interpolation method. It uses the polynomial of degree k—1, inter- 
polating k successive ordinates, as an approximation to F(x) only within a unit 
interval in x, centrally located with respect to its k defining ordinates. Assuming k 
fixed, successive approximating arcs for F(x) are thus obtained which present dis- 
continuities on passing from one arc to the next if k is odd, or discontinuities in their 
first derivatives if k is even (see section 2.121). Actually these discontinuities are 
irrelevant in our present case of an analytic function F(x). Indeed, if the interpolated 
values obtained are sufficiently accurate, these discontinuities will be apparent only 
if we force the computation beyond the intrinsic accuracy of the y,. 

The situation is quite different if y, are empirical data. In this case we are to 
determine an approximation F(x) which, for x=, may disagree with y, by amounts 
depending on the accuracy of the data, provided we thereby improve the smoothness 
of the resulting approximation F(x). In various applied fields such as Ballistics and 
Actuarial mathematics it is at times desirable to compute very smooth approxima- 
tions F(x) to an accuracy surpassing by far the accuracy to which the physical or 
statistical function involved may be determined. This physically unjustified accuracy 
becomes desirable whenever the approximation F(x) enters into numerical processes 
of some complexity, such as the numerical solution of differential equations. Modern 
electronic computing machines, especially, require a good amount of forced mathe- 
matical accuracy in such auxiliary tables in order to avoid the excessive accumulation 
of rounding errors in the computation of the solution. These remarks justify the de- 
sirability of approximation methods to empirical data furnishing easily computed 
approximations F(x) which are very smooth functions of x. Approximations meeting 
these requirements are of two kinds: 1. Polynomial approximation, where F(x) is com- 
posed of a succession of polynomial arcs meeting with a certain number of continuous 
derivatives. 2. Analytic approximations, where F(x) is an analytic and regular func- 
tion of x for all real values of x. 


* Received Oct. 18, 1945. 
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Important work concerning polynomial approximations is to be found in the ac- 
tuarial literature under the subject of osculatory interpolation. Of the extensive litera- 
ture we mention especially the fundamental work of W. A. Jenkins and the valuable 
systematization of the subject by T. N. E. Greville.! Especially important are those 
formulae derived by these authors which do not strictly interpolate the given ordi- 
nates, but rather combine the operation of smoothing the data and the operation of 
interpolation in one formula. Mr. Jenkins discusses interpolation formulae written 
in the convenient Everett (or Steffensen) form. Mr. Greville’s starting point is his 
elegant expression of each polynomial arc in terms of the end point values of those 
derivatives which are to be continuous on passing from one arc to the next. Each of 
these two modes of attack has its peculiar advantages and one or the other seem 
indispensable for an algebraic treatment of the subject. The present writer has found 
the Lagrange form (explicitly in terms of the ordinates y,) of such formulae preferable 
for two reasons: 1. The Lagrange form seems better adapted to computation with 
modern desk computing machines and undoubtedly superior for computation with 
punch-card machines. 2. The Lagrange form suggests a treatment of the subject by 
means of elementary concepts of Fourier analysis which, firstly, affords a more ex- 
haustive treatment of the problem of polynomial approximations, secondly, shows 
how to extend these methods so as to furnish analytic approximations. 

The explicit Lagrange form of the k-point central interpolation method, as well 
as of all the interpolation formulae of osculatory interpolation, is extremely simple in 
its formal appearance. Indeed, to every such formula corresponds an even function 
L(x), defined for all real values of x, in terms of which the corresponding formula may 
be written as follows 


F(x) = by Ynl(x — n). (1) 
The simplicity of this formula springs from the fact that it depends on the single func- 
tion L(x) which describes the formula completely. Incidentally F(x) =L(x) if 


y= 1, m=O (n¥0). (2) 


Thus every interpolation method of this kind exhibits its corresponding L(x) if we 
apply the method to the ordinates (2) (for an example see section 2.121). 

The polynomial interpolation formulae arise from (1) if L(x) is a composite poly- 
nomial function of arcs defined by various polynomials in successive unit intervals, 
such that L(x) =0 for sufficiently large values of | x| (for an important example see 
chapter II, formula (11)). The number of continuous derivatives of F(x) is, of course, 
equal to the number of continuous derivatives of L(x) for all real x. 

We obtain the formally simplest interpolation formula (1) if we choose 

sin rx 
L(x) = ——— (3) 
WX 

1W. A. Jenkins, Osculatory interpolation: New derivation and formulae, Record of the American In- 
stitute of Actuaries, 15, 87 (1926). 

Thomas N. E. Greville, The general theory of osculatory interpolation, Transactions of the Actuarial 
Society of America, 45, 202-265 (1944). 

W. A. Jenkins wrote four papers on this subject of which the above paper is the first. References to 
the other three papers are found in the excellent bibliography in Greville’s paper. 
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in which case (1) becomes 
: sin r(x — n) 
F(x) = pi a Se (4) 


Nn=—0o a(x _ n) 


This expression which interpolates the ordinates y,, is known to mathematicians un- 
der the name of the cardinal series.* For this reason we wish to call the general formula 
(1) a formula of the cardinal type, referring to L(x) as the basic function of the 
formula. 

The aim of the paper, of which the present Part A is the first, is twofold. Firstly, 
we propose to carry through to a certain stage of completion the important actuarial 
work concerning polynomial approximations. Incidentally, our work will answer Mr. 
Greville’s conjecture (loc. cit. pp. 212-213) concerning the existence of an “ordinary” 
interpolation formula furnishing an approximation F(x) composed of polynomial 
arcs of degree m+2, having m continuous derivatives and such that if the data y, 
are the values of a polynomial of degree m—1 then F(x) reduces identically to that 
polynomial. In Part B it will be shown how to obtain such formulae for every value 
of m. (The case of m=2 reduces to Jenkins’ formula mentioned in section 2. 122.) 
Secondly, we shall derive formulae of the cardinal type (1) with basic functions L(x) 
which are analytic and regular for all real or complex values of x. The classical basic 
function (3) is of course analytic; however, its excessively slow rate of damping, for 
increasing x, makes the classical cardinal series (4) inadequate for numerical pur- 
poses. Our analytic L(x), derived in chapter IV, dampen out exponentially. In Part B 
we will derive similar L(x) which will dampen out even faster: like exp(—C?x*). 

The paper is divided into five chapters. In chapter I we discuss the general prob- 
lem of smoothing by means of a linear compound formula. This discussion, by no 
means exhaustive, is to serve as a guide to what is likely to be useful among formulae 
of the cardinal type (1) which smooth and interpolate at the same time. It serves to 
restrict somewhat the arbitrariness of the problem. The rather obvious idea of us- 
ing cosine polynomials (or series) in this connection affords the possibility of a brief 
exposition of this subject in the more scientific manner of E. de Forest, W. F. 
Sheppard, E. T. Whittaker, and others, and may be followed up elsewhere. 

Chapters II and III form the common foundation of both parts A and B. In chap- 
ter II we describe the interpolatory properties of the formula (1) in terms of extremely 
simple properties of the Fourier-transform 


g(u) -{ L(x) cos uxdx (5) 


of the basic function L(x) (Theorem 4). Thus we are assured that our formula 
(1) will be exact for (i.e., reproduce) polynomials of degree k—1, provided g(u) —1 
has a zero of order k at u=0 and g(u) has zeros of order k at all points u=27n 
(n=+1, +2,---). This elementary fact is reminiscent of N. Wiener’s fundamental 

2See J. M. Whittaker, Interpolatory function theory, Cambridge Tracts in Mathematics, 1935, pp. 
62-64, for a discussion of the relation between the cardinal series and Stirling’s interpolation series. The 
cardinal series was probably first investigated in an important mémoire by Ch. J. de la Vallée Poussin, 
Sur la convergence des formules d'interpolation entre ordonnées équidistantes, Bull. Acad. Roy. Belgique, 
1908, 319-410. 
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description of the closure properties of the family of translation functions {L(x«—))} 
in terms of the zeros of g(u). Chapter III contains a somewhat general discussion of 
polygonal lines, the individual arcs of which are polynomials of degree k —1, joined to- 
gether with k—2 continuous derivatives. A general parametric representation of such 
curves is obtained (Theorem 5) which greatly facilitates their use for the purpose of 
approximation of data. For k=4 they represent approximately the curves drawn by 
means of a spline and for this reason we propose to call them spline curves of order k. 
These polynomial spline curves are finally smoothed out, by means of one-dimen- 
sional heat flow during the time interval ¢t, into analytic spline curves of order k. An 
analytic spline curve of order k is represented by a series of the cardinal type 


x 


F(x) = >> faMi(x — n, 2), (6) 
where the basic function M;,(x, ¢) is defined as 


_s* ,/2 sin u/2\* 
M,(x, t) = —{ e~ tu) (- —]} cos uxdx, (7) 


4 u 





while the coefficients f, may be thought of as arbitrary parameters. 

The family of functions (6) forms the basis of our work. Its principal advantages 
for purposes of numerical approximation spring from two sources: 1) The basic func- 
tion M;,(x, t) dampens out like exp(—<x*t~') (see III, formula (39)). As seen from our 
Table I, for k=4 and t=0.5, we have M,(x, 1/2) =0 to something like 10 decimal 
places for |x| 25. This causes the great flexibility of the graph of F(x) on varying 
the parameters f, and the ease in computing F(x). 2) The family (6) contains (or 
represents) all polynomials of degree k—1. The simplest analytic family of this type 
is obtained for k=Q and ¢>0 when (6) becomes 


F(x) = DU fn sig gran, (8) 
n=—co V mt 
This family obviously still enjoys the first property. However, (8) fails badly in its 
ability of representing even the simplest types of curves because of the low value of 
k=0. Indeed F(x) =0, for all f,=0, is the only constant value (8) is capable of repre- 
senting. 

Chapter IV contains the chief results of the present Part A. We show how the 
family of curves (6) can be used to approximate given data. First we derive an ana- 
lytic interpolation formula of the cardinal type (1) which leaves the given ordinates 
unchanged (Theorem 8). Secondly we extend the result to a family of formulae de- 
pending on a positive smoothing parameter e such as to combine a certain variable 
amount of smoothing (depending on e) with the operation of interpolation (Theorem 9). 

In collaboration with Lt. J. H. Levin, the author has had the opportunity of ap- 
plying on a large scale this analytic approximation method at the Ballistic Research 
Laboratory, Aberdeen Proving Ground, Maryland. The computations were per- 
formed on punch-card machines. The given equidistant data y, were the values of 
the drag coefficient of a projectile as a function of its velocity. Since very accurately 
computed values of the derivatives F’(x), F’’(x) of the approximation F(x) were also 
desired, it seems doubtful if any of the existing osculatory interpolation formulae 
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would have furnished satisfactory results in view of the complicated trend of the data 
to be approximated. 

In the last chapter we discuss procedures for the accurate computation of the func- 
tions and constants tabulated at the end of the paper. The most noteworthy problem 
encountered in this connection is the following: Let 


F(z) = >> ans” (9) 


be a Laurent series which converges ina ring a<|z| <8. We assume furthermore that 


f(z) does not vanish in this ring: 
F(z) #0, (a<|z| <8). (10) 


Under these circumstances we have an expansion of the reciprocal 





= > wns”. (11) 


If the coefficients a, of the expansion (9) are given numerically the problem consists 
in finding very accurate numerical values of the w,.’ A very efficient iteration method 
solving this problem has been developed by H. A. Rademacher and the author. It 
solves the similar problem of finding the expansion of the mth root of f(z) and gen- 
erally of any algebraic function of Laurent series. This subject will be discussed else- 
where in a joint publication with Professor Rademacher. 

In a sequel to these papers we expect to discuss the fitting of curves of the form 
(6) to data, in the sense of least squares. This will be accomplished by constructing 
series of the cardinal type (1) which also enjoy the orthogonality property 


ss 1 if x=0 
f L(x)L(x — n)dx = { ; 
0 if n#0. 


—o 


This construction reduces to the problem of computing the Laurent expansion of the 
square root 4/f(z) of an expansion (9). 

The author wishes to express his appreciation for the encouraging interest shown 
in his work by Dr. A. N. Lowan of the Mathematical Tables Project. He has bene- 
fited much by the helpful advice of Dr. L. S. Dederick, Major A. A. Bennett, Lt. J. H. 
Levin and others. Especially valuable were the author’s frequent discussions with 
Dr. C. B. Morrey. The tables were computed by Mrs. Mildred Young. The author 
takes this opportunity of expressing his thanks to the officials of the Ballistic Re- 
search Laboratory for their permission to publish these tables. 

The reader who is mainly interested in the numerical applications, may pass di- 
rectly from this point to the Appendix, where the use of the tables is fully explained 


and one example is worked out. 


8 It should be remarked here that also A. C. Aitken’s computation in 1925 of the coefficients of E. T. 
Whittaker’s smoothing method amounted to the expansion in a Laurent series of a certain simple rational 
function. See E. T. Whittaker and G. Robinson, Calculus of observations, London and Glasgow, 1940, 
pp. 308-312. 
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I, DEFINITIONS OF SMOOTHING AND SMOOTHING FORMULAE 
1.1. A definition of smoothing formulae. Let {y,} (n= --- —2, —1,0,1,2,---) 
be a given sequence or “table” which we wish to smooth. This smoothing operation 
is ordinarily performed by means of a formula of the following type 


Fe = Fa~ylep + “.% + Vn—1L1 + VnLo + Yn+1b-1 + = + YatpL—p, (1) 


where the numerical coefficients L, are symmetric about the middle term Lo, i.e., 
L,=L_,. The linear transformation (1) if applied to the original sequence {yn} will 
transform it into the smoothed sequence { F,}. By extending the definition of L,=0 
for | »| > p we may rewrite (1) as 


If y, =const.=c, we also wish that F,,=c; therefore 
pe a (3) 


is a natural requirement. 

When does the formula (1) actually smooth? As an example let p=1 and let the 
coefficients L, be (—1, 3, —1). If we now apply the formula (1) to the periodic se- 
quence 
1ynf = | --+,0,1,0,1,0,1,--- 
we obtain 

{F,} = {---,—2,3,-—2,3,-2,3,--+} 


which is a good deal rougher than the original sequence. Obviously this situation de- 
serves some clarification. 

There seems no doubt that the “smoothness” of a sequence {y,} depends in some 
way on its differences of higher order, especially on the sums of their square. We also 
notice that the formula (2) agrees with the rule of multiplication of Fourier series. 
This suggests the use of such series. 

Let us assume for the moment that 


De | yn| < @. (4) 
We now define a function J(u) by 
T(#) = 2: aie (5) 


r= oc 
and call it the characteristic function of our sequence {y,}; it is a complex-valued con- 
tinuous function of u of period 27. 


Now (5) implies 
e~*T(u) = Z. yreiir—Du = DS) yay seinu 


and by subtracting (5) we get 


(e-* — 1)T(u) = > Ave. (6) 


t=—s 


This shows that we obtain the characteristic function of the sequence {Ay,} of first 
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differences of {y,} by multiplying the characteristic function T(u) of {yn} by the 
factor e~*“—1. Generally 


- 4 


(e-iu — 1)"7T(u) ~ p> A™y,e'"" (m = 0, 8 x & oes + (7) 
Since |e~'“—1| =2 |sin(u/2)|, the Parseval relation furnishes the equation 
20 1 Qn ’ 
> (A™y,)? = " f [2 sin (u/2)]?™| T(u) [2du (m = 0). (8) 
n 0 PA 0 


These formulae furnish an expression for the sums of the squares of the differences 
of any order in terms of the characteristic function 7(u) of the sequence. 
Let us now turn to the “smoothed” sequence { F,,}. Let 


o(u) = } L,e*™" = Lp + 2L; cos u + 222. cos 2u+:-- (9) 
be the characterstic function of the sequence {L,}. We shall also refer to ¢(u) as 
the characteristic function of the smoothing formula (2). Notice that ¢(u) is always real 
and even. By multiplication of the two Fourier series (5) and (9) we obtain, in view 
of (2), 

T(u)o(u) = >> Fyein. (10) 

n=—c 

Hence the characteristic function of the “smoothed” sequence { F,} is obtained by 
multiplying the characteristic function T(u) of {y,} by the characteristic function 
o(u) of the smoothing formula (2). By now applying (8) to the sequence { Fr} we 
obtain 


- 9 


1 pe see 
D (a"r,)? = - f (2 sin (u/2))2™| T(1) |2(b(u))2du, (m=O). (11) 
‘ 2rd o 


A comparison of the relations (8) and (11) will readily furnish an answer to the ques- 
tion: what is a smoothing formula? Indeed, we notice that the integrands in (8) and 
(11) differ only, for each fixed value of m, by the factor @(u)* in (11). This justifies 


the following definition. 


DEFINITION 1. Let L, be a symmetric sequence of coefficients, t.e., L_»=Ln. The 


formation of the weighted means 


F,= >> vlna» (n =0, +1, +2,---) (12) 


is said to be a smoothing formula if 


>, L. = 1, (13) 
L|L.| < @, (14) 


while the characteristic function 


o(u) = >> Lyei™ = Lo + 2L, cos u + 2L2 cos 2u+--: (15) 


—*” 
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satisfies the condition 
—1< ¢(u) <1, (OS u S 2n). (16) 


The necessity of the condition (16) is justified as follows: By a comparison of (8) 
and (11), in view of (16), we obtain the inequalities 
2 o) 


z (A"F,)? < be (A™y,)?2, (m = 0,1, 2,--- ). 


n=—o n=—2o 


Actually the equality sign in one of these relations will arise only under highly ex- 
ceptional or else trivial conditions. This remark should make it clear why the smooth- 
ing quality of a formula violating (16) should be highly questionable. 

So far we were concerned merely with the ability of a formula (2) to smooth the 
sequence. However, the discrepancies between the two sequences also deserve atten- 
tion. By subtracting (10) from (5) we obtain 


T(u)(1 — o(u)) = =. (Vn on F,,)e*™ 
and therefore 


>» (vy, — F,)? = -f | T(u) |?(1 — o(u))2du. (17) 
0 


2r 


n=—2O 


A comparison of the integrands of (17) and (11) reveals the obvious fact that strong 
smoothing may be achieved only if we allow relatively large discrepancies between 
F,, and y,. Indeed, the integral of (17) will be small only if @(u) differs but little from 1, 
while strong smoothing requires as small a $(u) as possible. 

1.11. Examples of smoothing formulae. (a) Our trivial example Lop =3, Zi: =L_1=—1, 
L,=0 (n>1) has the characteristic function ¢(u) =3—2 cos u. We find ¢(u) 21, with 
(7) =5, which rules it out as a smoothing formula. 

(b) If L,20 for all n, and }OL, =1, then (12) is always a smoothing formula. 
Indeed 

| ¢(u)| = > Lei" | < 2. L,| = 1. 
Thus 
Fes +50 + edi (18) 


is a smoothing formula with 
o(u) = (1+ 2 cos u)/3. 


Let 


o;(u) = | (1) | =|{1+ 2 cos u | /3= Oe gy COs nu. 


Since (¢(u))?=(¢1(u))? it is clear from (11) that the formula (18) and the formula of 
characteristic function ¢:(u) have identical smoothing powers. However, since 
0<1—¢:(u) <1—¢(u) for 27/3 <u <47/3, we see by (17) that the formula 


,(1) (1) 
F, = > vba» 


will alter the sequence {yn} much less than (18) will. 
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(c) Generally, our formula (17) shows that it is desirable for an efficient smooth- 
ing formula to have its characteristic function satisfy the more restrictive condition 


0 < ¢(u) <1. (19) 


1.12. A comparison of smoothing formuiae. Again our relations (11) justify the 
following definition. 


DEFINITION 2. Let ¢:(u) and de(u) be the characteristic functions of two smoothing 
formulae. We say that the first is stronger than the second if 


| oi(u) | S| po(u) |. (20) 
with the inequality sign holding for some value of u. 


Later in this paper we shall set up a basic sequence of smoothing formulae of 
progressively greater strength according to this definition. Here we remark only that 
two smoothing formulae cannot in general be compared on the basis of this definition. 
However, the following remark seems obvious. Let 

F, = > yebn-» (21) 
” 
be a smoothing formula of characteristic function @(u). The iteration, or repetition, 
of (21) may be thought of as another smoothing formula and its characteristic func- 
tion is found to be (¢(u))?. Since | b(u)| <1 obviously 


(p(u))? < | o(u) |. 


This shows that the formula (21) and the sequence of its successive iterates form a 
sequence of smoothing formula of progressively increasing strength. 

1.13. Smoothing formulae which are exact for polynomial values of a given degree. 
The following definition is in common use. 


DEFINITION 3. A smoothing formula (2) is said to be exact for the degree m if it re- 
produces exactly the values {y,} of a polynomial of degree not exceeding m. 


If 


F, = 2) vole» (22) 
is to be exact for the degree m, it is obviously sufficient that it be exact for the basic 
monomials 1, x, ---,a«™. Thus the exactness for the degree m is equivalent to the 
relations 

- 
ilies ) vba» (s = 0, I,s°*, m). (23) 
v=—00 


Let us now assume for simplicity that the sequence of coefficients L, tends to zero 
exponentially as n— ©, i.e., we assume the existence of two positive constants A and 
B such that 





| La| & AeFim 


for all values of n. This implies that the function 
(u) = 7. L,e*™ 


is regular in the strip 
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| Iu| < B 
of the complex u-plane. Now 
o(u) = > Lie** = D> Lyne? e7ine 


and 


e*™“d(u) = r» e”*L,_». 
y 


We now expand both sides in ascending powers of u and compare like coefficients. 


Since 
2 4 


u u 
g(u) =1+— 9" +— 9M 4---, 
2! 4! 


we get the identities in n 


s S “= 
P = ( Jorn + (| ) oon ‘$_...= DO vl, (s = 0,1,2,---). 


J 
y=—oo 


A comparison with (23) will show that a smoothing formula is always exact for a 
highest degree which is always odd. It also proves the following proposition which 
may evidently be established under conditions less stringent than the ones we used. 


THEOREM 1. A smoothing formula (22) is exact for a degree 2v+-1 if and only if 
o(u)—1 has at u=0 a zero of order 2v+2, 1.¢., 


¢(0) = (0) = --- = 90) = 0. (24) 
As an illustration we mention the formula 
F, = = (— Yn-s + On-1 + 16yn + Int — Yuta) = Yn — - ity, — = 5°. (25) 
32 16 32 


of characteristic function 

o(u) = (8 + 9 cos u — cos 3u)/16. (26) 
We find that $’’(0)=0, hence (25) is exact for cubics. The symmetry property 
o(u) +(x —u) =1 shows that 

o(r) = ¢'(x) = $' (4) = (x) = 0. 
This results in rather strong smoothing power. The formula (25) is part of a sequence 
of formulae, the next one of this kind being 


1 
3 Va—5 25 yn—-s + 150yn-1 + 2560 + 150yn41 — 25 yn+3 + 3 Yn+5) (27) 


F, = — (33 
512 


or 


5 15 3 
F,, _ Yn + vphake: é°y,, + ee 5*y,, + — §}y,,. 
32 256 512 
Its characteristic function 
o(u) = (128 + 150 cos u — 25 cos 3u + 3 cos 5u)/256 


again enjoys the symmetry property $(u)+¢(7—u) =1. Also ¢(u)—1 has a zero of 
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order 6 at u~=0, hence (25) is exact for quintics, while ¢(u) has a zero of order 6 at 
u=m resulting in strong smoothing power. 

1.2. Smoothing a finite table. In 1.1 we have discussed the smoothing of an infinite 
table {y,} which is such that the series of the absolute values of its entries converges. 
By (8), (11) and the inequalities (16) we have found that the sum of the squares of 
the differences of order m is diminished by smoothing. This is true form=0,1,2,---. 
Now we shall discuss briefly the practically most important case of a given finite table 


{yn} (n = 0,1,---, N). (28) 


To fix the ideas we assume the following simplest concrete situation: the third 
differences A*y, are slowly varying and of slowly varying signs, while the A‘y, are of 
random signs. In this situation we naturally wish to minimize the 4th differences of 
the table. Now we form an average value of the A*y, at each of the two ends of the 
table and we extend the column of A*y, with the corresponding constant average value 
at each end.‘ Thus the A*y, are defined for all m having one constant value for n> N—3 
and another constant value for <0. Now we extend the definition of y, for all m from 
the values of the third differences. Also, we compute the A‘y, for the extended infinite 
table. Clearly A‘ty, =0 for n< —1 or n>N-—3. Let 


T,(u) = >> Atynein™ 


be the characteristic function of the sequence of 4th differences, the series containing 
really a finite sum of terms only. 
Let us now apply to the extended table y, a smoothing formula 

F, = Do volun (29) 
of characteristic function ¢(u), which is exact for cubics. The result is the new sequence 
1 Fa} (— © <n< ). Evidently y, are the values of cubics for large n| and therefore 
F,=y, for large |n|, hence also A*F, =A*y, and A‘F,,=0 for large ||. Notice also 
that we may think of the sequence {A‘F,} as arising from {A‘y,} by the smoothing 
operation (29). Therefore 





> (A*y,)? = —fi| T 4(u) |2du 


n=—oO 


and 





*6(u)*du. 


a 1 2r 
; 2 (A4F,,)? = —{ | T ,(u) 
T 0 


n=— 


Generally 








20 1 Qe 

> (Attmy,)? = —{ (2 sin «/2)?™| T4(u) |2du, (m = 0) (30) 
n 00 ~ 0 

Bal 1 2nr 

> (Att™F,)? = —f (2 sin u/2)?™| T4(u) |2o(u)2du, (m = 0). (31) 
n=—0 TH 0 


‘ Compare G. J. Lidstone, Note on the computation of terminal values in graduation by Jenkins’ modified 
osculatory formula, Transactions of the Faculty of Actuaries (Scotland), 12, 277 (1930). 
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A comparison of (30) and (31) shows that the sums of the squares of the fourths and 
subsequent differences have been decreased by the smoothing operation. No such 
statement can or should be inferred concerning the finite sums of relevant differences 
of orders 0, 1, 2, and 3. 


II. INTERPOLATION FORMULAE 
2.1. Interpolation formulae of the cardinal type. Let 
F, = 20 Ln (1) 


be a smoothing formula. If we apply it to the “elementary” table 


yo = 1, in =0 (n #0), (2) 


then F,=L,. The even sequence {Ln} may therefore be regarded as the smoothed 
version by (1) of the elementary table. Now suppose that we are given not only the 
even sequence of ordinates L, but an even function L(x) defined for all real x and such 
that L(n)=L,. Then we may replace the integral variable in (1) by the continuous 
variable x and we obtain the formula 


00 


F(x) = Zz. y L(x — v). (3) 

We call L(x) the basic function of the formula (3). The chief aim of this paper is 

to point out that the subject of interpolation is truly dominated by the formulae of 

the type (3), the kind of approximation desired depending only on the choice of the 
basic function L(x). The particular basic function 


sin rx 








nx 
gives rise to the series 
- sin r(x — pv) - 
F(x) = V,——— (5) 
Rea a(x — v) 


which is well known to mathematicians and referred to as the cardinal series. For this 
reason we wish to call (3) a series, or formula of cardinal type. 

We notice here for further reference that the basic function (4) may also be writ- 
ten as a Fourier integral as follows. 


r 


sin 7x 1 


——————“(Ss— e**"du. (6) 

1x 2nd u,; 

2.11. The two kinds of interpolation formulae, ordinary or smoothing. For integral 
values of x= our formula (3) becomes 


io) 


F(n) = > y,L(n — v). (7) 


v=—oO 


Equation (3) is an interpolation formula in the usual sense if F(n) =y,, for all m, and 
this is the case if and only if L(x) satisfies the conditions 
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L(0) = 1, L(n) =0 (n #0). (8) 


Otherwise, (7) isa smoothing formula. We shall follow the accepted actuarial practice 
of referring to (3) as an ordinary interpolation formula if (3) reproduces exactly the given 
ordinates {yn}. Otherwise we call (3) a smoothing interpolation formula. 

2.12. Examples of interpolation formulae of the cardinal type. Later in this paper 
we shall discuss various classes of such interpolation formulas all arising from a com- 
mon general theory. For purposes of orientation and illustration we mention here a 
few concrete examples. 

2.121. The k-point central interpolation formula. Let k be a fixed integer 
(=1,2,3,-- +). By k-point central interpolation we mean the interpolation method 
whereby the polynomial of degree at most k —1, defined by k consecutive ordinates y,, 
is utilized within an interval of unit length centrally located with respect to the set of 
defining ordinates y,. This set of k defining ordinates y, is shifted up by one unit in 
the subscript for interpolation in the next unit interval. It seems obvious that this 
kind of interpolation is performed for any real value of x by a formula of the cardinal 
type 


F(x) = > nC x(x — n) (9) 


n=—o 


with a suitable function C;(x). To obtain this function, it is sufficient to interpolate 
the elementary table (2) by means of this method of k-point central interpolation. 
The graphs indicate the resulting C;,(x) for k=1, 2, 3, and 4.5 It is found that C;(x) 


il 


C, (x) Co (x) 








1 
' 
H | 
4 ‘ 
‘ ' 
‘7 A 
yw 1 2 <2 o} 0 


—— 


“2 -1 -l2 


A V4 
SI. C3(x) 


“8 
! 
' 
; ' \ n 
ee 3/2 _. 





k 

a os. 
ol 1\ wae 
‘ 





5 These graphs indicate geometrically the construction of the successive arcs of these curves. Thus 
C;(x) is defined in the interval 1/2<x<3/2 by the parabola passing through the points (0, 1), (1, 0), 
(2, 0). Similarly C,(x) is defined in —1<x<0 by the cubic which takes the values 0, 0, 1,0atx=—2, —1, 
0, 1 respectively. We mention incidentally the following general analytic expression of the basic function 
C(x) of k-point central interpolation. In terms of the “central” factorial 


x(x? — 12)(4? — 22) --+ (x? — (» — 1)*) if k = 2p, 
xk)“ = 


(eB) Be RGB) tee 


we define the corresponding truncated function 
[k]-1 “a if s>0 
0 if «<0 (k = 1,2,3,-++). 


+ = 
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is a polygonal line composed of arcs of degree k—1. Also C;(x) itself is continuous with 
a discontinuous first derivative if k is even. For an odd k, C;(x) itself is discontinuous, 
the value assigned at a point of discontinuity being the arithmetic mean of the two 
local limits. Evidently for k=2 our formula (9) is identical with the method of linear 
interpolation and the graph of F(x), as given by (9), is identical with the polygonal 
line of vertices (”, yn). 

For further reference we mention the following formulae which are valid for all 


1 * sinu/2 ~ 
Ci(x) = — —— e'“"du, 
=< 


2r u/2 


1 ¢°/sinu/2\? | 
C2(x) = —{ ms Fee, 
2rJ_.\ u/2 
1 * /sin u/2\8 1 
C3(x) = —{ ( —— a aS uz edu, 
2B u/2 8 
1 * /sin u/2\4 1 
Cz) = f ee ) (1 ee ’) ei“zdu. (10) 
2arJ_.\ u/2 6 


The k-point central interpolation method is the most important method for the 
interpolation and the construction of tables of analytic and regular functions. How- 
ever, for the construction of tables of empirical functions, the low order of continuity 
of C;(x) is at times a serious limitation of this method. It seems indeed evident that 
the continuity properties of the linear compound 


F(x) = >> y,L(x — n) 


n 


real values of x 


are directly dependent on the continuity properties of the basic function L(x). We 
turn now to an interesting example of an “osculatory” interpolation formula having a 
basic L(x) enjoying stronger continuity properties. 

2.122. An osculatory interpolation formula of W. A. Jenkins. We define a basic 
function L(x) for x $0 by 


0 if xs —3 
1 . 
| —- — (2+ 3X2 + 2) if -—35x%5-2 
12 
L(x) = 4 ? (11) 
| gr et ne + Oe + 9 if —2sx*xs-1 
- 
< (x + 1)(6 — 6x — 9x? — x) if —1s x50 


oval temas: 


The definition of this function is to be completed for x=0 by continuity, if k is 
metic mean of the two limits, if k is odd. Then 


given, and by the arith- 


k [kj)-1 


+ (-» <a4< »), 


Cx) = ——— 
(k — 1)! 
where 6* is the symbol of the &th central difference of step unity (compare Theorem 3 of section 3.11). 
We will return to this subject in Part B where the extremely simple law of formation of the integrals (10) 


will also be given. 
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and extend its definition by the requirement L(—x)=L/(x) to all real values of x. 
Since the conditions (8) are visibly verified we see that 


F(x) = > ¥nL(x — n) (11’) 


is an ordinary interpolation formula. A closer inspection of the composite polynomial 
function (11) will show that L(x), L’(x), and L’’(x) are all continuous for all real val- 
ues of x. Using a customary mathematical terminology we may say that L(x) is of 
class C’’. Moreover, in various ways it may be shown that the formula (11’) is exact 
if the y, are the ordinates of a polynomial of degree 3 or less, i.e., F(x) becomes iden- 
tical with that cubic polynomial. 

It is of interest to compare Jenkins’ formula (11’) with the 4-point central inter- 
polation formula (9) (for k=4). Both are exact for cubics. C,(x) is of class C, while 
the present L(x) of class C’’. This was achieved by increasing the complexity of the 
basic function in two ways: 1) The interval where L(x) is non-vanishing was increased 
from |x| <2 to | x| <3. 2) The degree of the polynomial arcs has increased from 3 
to 4. Later Jenkins’ formula (11’) will appear as a member of a sequence of interpola- 
tion formulae of similar characteristics. Here we mention that the basic function (11) 


may be expressed in the form 


i * /sin u/2\4 sinu\ | 
L(x) = —{ (= *) 4+ cos u — 4— ])e‘*du (11’’) 
2rJ_.\ u/2 u 


for all real values of x. 
2.123. A smoothing interpolation formula of W. A. Jenkins. We define a basic func- 


tion L(x) by 





0 if «xs — 3, 
1 

— — (x + 3)3 if —3s x85 -2, 
36 

L(x)=4 1 (12) 
go re ee Bee? if —2s*s-1, 
Bis 
— (15 — 27x? — 1423) if —1isx<<0, 
18 
L(-— x) = L(x). 


This particular L(x), composed of cubic arcs, is of class C’. The formula® 


6 W. A. Jenkins writes his interpolation formula (11’) in the following Everett form 


e(g — 1 ME —1 
Fin + x) Vnt +- 5*y, a —_ iy So 





’ (Osxs1,x+é=1). 





+ Yayi% + 57 yn 41 

Likewise his formula (12’) takes the form 
_ 1) 3 
F(n + x) = yn& + 8%yn ao ~ % 5 


8 3 
+ Vn4irt +- 57 yn +1 «tliat J _ 5 Yn51 a - 
a 6 “36 
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F(x) = >> y.L(x — ») (12’) 
, 
corresponding to (12) is exact for polynomials of degree 3 or less. However, while 
(11) was an ordinary interpolation foimula, the formula (12’) is a smoothing inter- 
polation formula. Since 


L(0) = 15/18,  L(i) =2/18, L(2) = — 1/36, 


while L(n) =0 for n23, we see that for x= (12’) reduces to a smoothing formula of 
characteristic function 
o(u) = bs (15 + 4 cos u — cos 2u). 
18 

We readily verify that $’’(0)=0 and ¢(7) =5/9 S<¢(u) $1. Hence (12’) reduces for 
integral x= to a smoothing formula, according to our Definition 1. On comparing 
Jenkins’ two formulae (11’) and (12’) we notice that they are both exact for cubics, 
giving rise to curves of class C’’. Since (12’) is only a smoothing interpolation formula 
while (11’) is an ordinary interpolation formula, it has been possible to lower the de- 
gree of L(x) from 4 to 3. We finally mention that the function (12) may be expressed as 


1 * /sin u/2\4 / 4 1 
L(x) = —{ (= — — — cos u ) e*“*du. (12’’) 
2a —« u/2 3 3 


2.2. A general theory of interpolation formulae of the cardinal type. In this sec- 
tion we shall discuss various characteristic properties of interpolation formulae of the 
cardinal type in terms of the Fourier transform of the basic function L(x). This dis- 
cussion will provide a sufficiently broad foundation for the subsequent development 
of specific formulae in the latter part of this paper. 

2.21. Characteristic properties of interpolation formulae. Some of the following defi- 
nitions have already occurred in the previous sections. For convenient reference we 
include them in our present enumeration of properties of an interpolation formula 

00 
F(x) = >> »,L(x — »). (13) 
rm —20 

a. We say that (13) is an ordinary interpolation formula if F(x) interpolates ex- 

actly the given ordinates yn, i.e., if 
L(0) = 1, Liv) =0 (v0). (14) 


b. We say that (13) is a smoothing interpolation formula if for x =n (13) turns into 
a smoothing formula 


F(n) = >> yL(n — »). (15) 


The term “smoothing formula” is meant, of course, in the sense of our Definition 1, 


section 1.1. 
c. We say that (13) is exact fur the degree k—1 if the relation 


P(x) = > P(n)L(x — n) 


n=—oO 
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is an identity for any polynomial P(x) of degree at most k—1. The last condition is 
in turn equivalent with the k identities 
x = > n(x — n) (v=0,1,---,k— 1) (16) 
out of which it can always be recovered by means of suitable linear combinations. 
d. We say that (13) preserves the degree k—1," if for any polynomial P(x) of degree 
v<k—1 we have an identity 


>> P(n)L(«x — n) = P(x) + (a polynomial of degree < v). (17) 
Notice that the leading term of P(x) is not altered by (13). Again in terms of the mo- 
nomials x” we may say: (13) preserves the degree k—1 whenever the k functions 


- a 


0,(x) = y x n’L(x — n), (vy =0,1,---,k—1), (18) 


Ls 


are polynomials of the form 
Ox) = aw + ayx'+---+a4,, (v=0,1,---,k— 1). (19) 


e. We say that (13) is of degree m and of class C*, if the basic function L(x) is a 
polygonal line of polynomial arcs of degree at most m joining in such a way as to 
result in a function L(x) having pw continuous derivatives. In the sequel, the junction 
points will always be either for integral values x=n or else for x=n+1/2. Conse- 
quently no “condensation” of discontinuities will result by the formation of the linear 
compound (13). Hence the interpolation curve F(x) will again be of degree m and of 
class C*. As examples we recall the formulae (11’) and (12’) of W. A. Jenkins, which 
are both of class C” and of degree 4 and 3, respectively. 

f. A formula (13) whose basic function L(x) is composed of polynomial arcs will 
also be referred to as a polynomial interpolation formula. We shall say that it has the 
span s if the even function L(x) vanishes identically for x >s/2, but not for x>s’ with 
0 <s’<s/2. Thus the k-point central interpolation formula (9) is of span k, while both 
formulae (11’), (12’) of Jenkins’ have the span 6. For obvious practical reasons it is 
desirable to work with polynomial formulae having as small a span as possible. 

g. We say that (13) is an anaiytic interpolation formula if the basic function L(x) 
is analytic and regular for all real x. The original cardinal series (5) is an example of 
this type. Obviously no analytic formula can possibly have a finite span. The role of 
the span is taken over by the rate of damping of L(x) as x increases. For obvious prac- 
tical reasons it is desirable to work with analytic L(x) damping out as fast as possible. 

2.22. The characteristic function of the basic function L(x). It was shown in chapter 

7 This property of interpolation formulae seems to have been neglected so far. It represents an im- 
portant weaker form of the condition of exactness for the degree k —1. Compare T. N. E. Greville, loc. cit. 
pp. 210-211, for our slight departure from the standard terminology. Jenkins speaks of a modified inter- 
polation formula in case the formula is not ordinary. The term “modified” seems natural in view of 
Jenkins’ construction of such formulae by modifying certain terms of Everett’s formula (the author is 
indebted for this last remark to Mr. Chalmers L. Weaver). It seems, however, less desirable if their con- 
struction is, as here, otherwise performed. 
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I that various properties of a smoothing formula 


Pie >, tha 
, 


are readily expressible in terms of its characteristic function 


¢(u) = > Le. 


Likewise, the properties of the interpolation formula 


F(x) = 7. y L(x — v) (20) 


? 


will largely depend on the behaviour of the function 


g(u) -{ L(x)e*“*dx -{ L(x) cos uxdx. (21) 
This even function g(u) is the Fourier transform of L(x). Following a terminology 
used in probability theory we shall refer to g(u) as the characteristic function of L(x). 

Under certain general assumptions which will always be verified in our applica- 
tion, the relation (21) may be inverted’ to 


1 es) 
L(z) = —{ g(uje*du, (22) 


Tf —« 


However, it should be remarked that at times our integrals are not absolutely con- 
vergent and that they then converge only as a principal value in the sense of Cauchy: 
lima ales An example of this kind is our first formula (10) 


1 °@sinu/2 - 
C(x) = — f 7 - ei@zdy, 
—. 





Changing u and x to 2ru and x/2m respectively we see 
x °sinwu - 
Ci{—] = — e'“tdy, 
2r -«o TU 
Inverting this relation we obtain 


sin rx 1 “s u\ . 1 
— = — Cy “) e“zdyu = — e“zdy (23) 
ye 2 


Wx ar 2a J _-» 








which is identical with (6) and shows that 


1 if < 
” ( if |ul<-r 
Pe c,(*) ‘i 
2r 


is the characteristic function of the original cardinal series (5). It is precisely the dis 
continuity of its characteristic function which causes the extremely slow damping of 


the basic function (23). (See 2.21, g.) 


1 if |ul=-r 


lo if |ul>-r 


8 See e.g. S. Bochner, Vorlesungen tiber Fouriersche Integrale, Leipzig, 1932, Satz 11b on p. 42. 
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Similar reasons of slow damping will rule out the following rather obvious method 
of turning a given smoothing formula 


F(x) = >> vba» 


into a smoothing interpolation formula 


F(x) = >> y,L(x — »). 
From 
o(u) = >> L,ein 
we derive 
1 r 
L,=— o(u)e*™du. 
\. 


Now we simply define a basic function L(x) by 


1 r 
L(x) = —{ o(uje*"du. 
2rJ_, 


The corresponding characteristic function g(u) is found to be 


o(u) if |ul<- 
g(u) = { 
if | u| > =. 


Again the discontinuities of g(u), or of one of its higher derivatives, will imply that 
the damping of L(x) is too slow for numerical purposes. Indeed, by partial integra- 
tions, L(x) is found to tend to zero as a certain negative power of x only, as x tends 
to infinity. (Concerning the order of magnitude of Fourier integrals for large values 
of x, see the theorem on page 11 of Bochner’s book quoted in our footnote 8.) 

2.23. Fundamental criteria in terms of characteristic functions. We shall now re- 
strict ourselves to basic functions L(x) which are everywhere continuous with the 
exception of possible “discontinuities of the first kind” (such as were exhibited by the 
basic functions L(x) of section 2.121). Moreover, we shall assume that L(x) dampens 
out exponentially. This means that we assume the existence of two positive constants 
A and B such that the inequality 


| L(x) | << Ae-Blzl (24) 


holds for all real values of x. This clearly rules out the basic function (23) of the cardi- 
nal series. The assumption (24) implies that the characteristic function 


g(u) = fi cenetas (25) 


is analytic and regular not only on the real u-axis but also in the infinite strip 
|Iu| <B (26) 


of the complex u-plane. It also implies that the expression 
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Pp 
lim Zz. g(u + 2nn)e?*in2 (27) 


2 n=—p 
converges uniformly in a circular neighborhood of «=0 and for every real value of x. 


The following theorem will demonstrate the usefulness of the characteristic func- 
tion of an interpolation formula. 


THEOREM 2. Let the basic function L(x) satisfy the condition (24). Let the correspond- 
ing interpolation formula be 
F(x) = >> L(x — »). (28) 


For integral x =n (28) reduces to the smoothing formula 
F(n) = >> y,L(n — »). (29) 


A. The characteristic function @(u), of the smoothing formula (29) is given by the 
relation 
o(u) = >> g(u + 2m). (30) 
In particular (28) is an ordinary interpolation formula (see 2.21, a) if and only if 
> g(u + 2nv) = 1. (31) 
B. The formula (28) is exact for the degree k—1 (see 2.21, c) tf the following two 
conditions hold simultaneously: 


g(u) — 1 has a zero of order k atu = 0, (32) 


g(u) has zeros of order k for all non-vanishing integral multiples of 2x;u=24n (n¥0). (33) 


C. The formula (28) preserves the degree k—1 (see 2.21, d) af the condition (33) holds, 
together with the additional condition 


g(0) = 1. (34) 


Remark. For some applications it is important to notice that an ordinary inter- 
polation formula which preserves the degree k—1 is automatically exact for the degree 
k—1. This seems evident a priori. It is also evident in terms of our criteria, for (31) 
implies 
g(u) —1= —> glu + 2rv) 


¥#0 


and the right-hand side has a zero of order k at u=0 by (33). 
Proof of A. Our formula (25) implies 


1 oo 1 (2p+1)r 
L(x) = — g(u)e**“du = lim —{ g(uje*“du 


20 Su poo le (2p+1)x 


12 fr 
lim — > g(u + 2mv)ei*e?**“=dy 


~-2 2n v=u—p —" 
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1 r Pp 
lim — { > glu + danetrind e'“*du, 


| de 2x —« \p=x—p 


1 T a) 
—{ { > g(u + dant} ei“zdy, 
2x J _s 


b==—cO 


II 


In particular, if «=n is an integer, we find 


1 T a) 
L(n) = —{ { es g(u+ ano) ened, (35) 
41 v/ _# —2 


Since the characteristic function @(u) of (29) is by definition the function of 
Fourier coefficients L(m) (see 1.1, (9)), the relation (30) is established. 
Proof of B. We wish to apply Poisson’s summation formula® 


> fx-—n) = > etrins f f(v)e-** "dv (36) 


to the function 


f(x) = e-**L (x). (37) 
By (37) and (25) we find 


f S (ve? "dy -f L(v)eiut2e—)edy = g(u + 2xn) 


hence by (36) 


Ro) oa 
en izu ) eT (x —n) = > g(u + 2en)e2*inz 


and finally 
pm e“ "L(x — n) = ett > g(u + 2rn)e****, (38) 


This identity actually holds for all real x and all real or complex values of u within 
the strip (26). It contains implicitly all the statements of Theorem 1. Thus forx=0 
it reduces to (30). To prove our statement B we assume x fixed and regard both mem- 
bers of (38) as functions of u, which we expand in series of powers of u, then equating 
the respective coefficients on both sides. On the left-hand side we have the expansion 


i) au” co) 


> a p n’L(x — n). 


yr Vi ne—co 


On the right-hand side, our assumption (33) implies that the terms g(u-+-2mn)e***"* 
(x0) do not contribute any terms in u of order less than k. Thus our identity (38) 


becomes 


r-) 7” » oo 
> — >> w(x — n) = e***g(u) + u* (regular function of x). (39) 
y=0 Vi- ne—oo 


On the other hand our assumption (32) amounts to 


» See S. Bochner, loc. cit., theorem 10 on page 35. 
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g(u) = 1 + w* (regular function). 


This and (39) imply 


© au’ ce) i) a” . 
= om 2. n’L(x —n) = >» a ap u* (regular function). (40) 
yp=0 Vi n=e—o | ea 


A comparison of the coefficients of the first k terms on each side of (40) furnishes the 
identities (16). This concludes a proof of B. 
Proof of C. Since g(u) is regular at u=0, and even, it has in view of (34) an ex- 


pansion of the form 
ae = a4 ay 
g(s) @ 1 — — *® + — ot — — ee t:::. 
2! 4! 6! 


We now define a sequence of polynomials by means of the generating function 


* (1u)” 
e*"9(4) = z=. Q,(x) - (41) 
v=0 V. 
or 
fe 4) u” 
e*"g(u/t) = Zz; Q,(x) - : (42) 
y= Vv: 
where 
’ a: , a F 
RIP See eee +t (43) 


A comparison of terms on both sides of (42), using (43), shows that 


QO,(x) = x’ + ( ; ) a2x”—* + (7 ) ove +---, (44) 


On substituting the expansion (41) into the right-hand side of (39) and by compari- 
son of the first k terms on both sides we find that the identities (18) and (19) are 
established. This completes the proof of our theorem. 

As a brief illustration of our criteria let us consider again Jenkins’ smoothing in- 
terpolation formula (12’) of 2.123. Its basic function is 


1 * /2 sin u/2\4 / 4 1 
L(x) = —{ (— :) (- — — — COs “) emrdu, (45) 
2nd _. u 3 3 


A simple method of evaluating explicitly such integrals in order to find the polynomial 
expressions (12) will be discussed later. An inspection of the characteristic function 


2 sin u/2\4 
g(u) = Cae) (4 — cos u)/3 
u 
immediately reveals that our condition (33) is verified for k=4. Direct expansion 
shows that g(u) =1—(7/240)u‘+ --- and (32) is also verified for k=4. The inter- 
polation formula (12’) is therefore exact for cubics. Also the fact that L(x) is of class 
C” is revealed by an inspection of the integral (45). Indeed we notice that g(u) van- 
ishes for u= © like u~‘. This implies that we may differentiate (45) twice under the in- 
tegral sign and that the integral 
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1 oo 
L(x) = —f g(u)(iu)e**du 
ond 
is also continuous since it converges absolutely. 


III. THE THEORY OF SPLINE CURVES 


The previous chapter provides a formal theory of interpolation formulae in terms 
of a basic function L(x) which, as yet, is largely arbitrary. The present chapter will 
furnish the foundation for the derivation of special basic functions which are readily 
computed with great accuracy and lead to interpolation formulae enjoying the proper- 
ties described in the previous chapter. 

3.1. Polynomial spline curves of order k. A spline is a simple mechanical device 
for drawing smooth curves. It is a slender flexible bar made of wood or some other 
elastic material. The spline is place on the sheet of graph paper and held in place at 
various points by means of certain heavy objects (called “dogs” or “rats”) such as 
to take the shape of the curve we wish to draw. Let us assume that the spline is so 
placed and supported as to take the shape of a curve which is nearly parallel to the 
x-axis. If we denote by y=y(x) the equation of this curve then we may neglect its 
small slope y’, whereby its curvature becomes 

1/R = y"/(1 + yy” 
The elementary theory of the beam will then show that thecurve y=y(x) is a polyg- 
onal line composed of cubic arcs which join continuously, with a continuous first 
and second derivative.'!® These junction points are precisely the points where the 
heavy supporting objects are placed. 

3.11. Description of spline curves of order k. Our last remark suggests the following 
definition. 


DEFINITION 4. A real function F(x) defined for all real x is called a spline curve of 
order k and denoted by I1,(x) tf it enjoys the following properties: 

1) It is compressed of polynomial arcs of degree at most k—1. 

2) It is of class C*~*, t.e., F(x) has k—2 continuous derivatives. 

3) The only possible function points of the various polynomial arcs are the integer 
points x =n if k is even, or else the points x =n-+1/2 if k is odd. 


Thus a II;(x) is a step function with possible discontinuities at the points 
>=n+1/2. A IIe(x) has an ordinary polygonal graph with vertices only at the in- 
teger points x=n. A II«(x) corresponds to the elementary mathematical description 
of an ordinary (infinite) spline with the “dogs” placed at all or only some of the points 
with x=n. 

It should be noticed that if a II,(x) is of class C*—', then II“- (x) must necessarily 
be constant for all x. Thus such a II;(x) reduces to a polynomial of degree k—1. It 
is just this relaxation of the requirement of the continuity of the (k —1)-order deriva- 
tive of II,(x) which turns the spline curve into a flexible and versatile instrument of 
approximation. Likewise, only the “dogs” (or “rats”) enable the ordinary spline to 
trace curves differing from the graph of a cubic polynomial. 

The special importance of spline curves will be due to the fact that by the addi- 


10 The author is indebted for this suggestion to Professor L. H. Thomas of Ohio State University. 
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tion of several spline curves of successive orders we may get any desired polygonal line 
of given degree m and class C“. 

3.12. The evaluation of certain Fourier integrals. Our further work is based on the 
consideration of the functions 


1 ° /2 sin u/2\* | 
M;(x) = — f ———} e*“*du (he 1,.2,05+3 = © Sae< oo), (1) 
vo on u 


They have been evaluated explicitly for low values of k by various authors." The fol- 
lowing general explicit representation is essentially due to Laplace (see J. V. Uspen- 
sky, Introduction to mathematical probability, 1937, Example 3, pp. 277-278). 
THEOREM 3. Let k be a positive integer. Define the function x*>' by 
k-1 r. if x20 
x4, = 
0 if «<0. 


+ 


For k=1 and x=0 this definition 1s modified to 0&1 =1/2. 
The following identity holds for all real values of x 


1 © 7/2 sin u/2\* 1 ee 
- -f ( ———— } e“"dy = ———-§ %, , (3) 
2r J _« u (k — 1)! 


where 5* stands for the usual symbol of the kth order central difference of step equal to 
unity. 


The identity (3) is correct for k=1. Indeed it is well known that 


1 © sin u 
~- du = 3. 
v0 8 


On replacing u by ux we get that 





: r 4 if z«>0 
1 * sin xu j ; 
—{ —— du = 0 if «=0 
dn J _« u : 
—} if «<0, 


or 


| = 


“ sin xu 0 
. ae a ae 
du = %4 A. 
24 J _« u 
Therefore 


1 * 2sin u/2 1 * sin (x + 3)u 1 * sin (x — $)u 
———— cos uxdu = — ee BE me ——__—__——— du 


2rd _» u Zid —« u Qa J _.n u 
1. f* sin ou cal 9 
= —6 — du = 6(x, — 4) = 6x,, 
2r me u 


11 See S. Bochner, Fourier analysis, Princeton University lectures, 1936-1937, where our M,(x) are 
worked out for k=1, 2, 3, and where the increasing smoothness qualities of these functions are clearly 
noted. Bochner also considers the integrals 

o3.. 1 ¢* /2 sin u/2\* sin ux 
Mi(x) = — — du 
2r u u 





Using Theorem 5 below, we readily obtain the identity M,*(x) = —}+(1/k!) ox. 
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and (3) is established for k=1. 
Let us consider for the moment the sequence of functions 
1 i bes 
Ni(x) = ——— 5b x, . (4) 
i-th " 
We have already shown that 
M(x) = N,(x) (5) 


holds for k=1. Assume now that (5) holds. We wish to show that the similar identity 
for k+1, rather than k, arises from (5) by performing the operation 


+ aie 
z—1/2 


on both sides of (5). This will be accomplished if we prove that 


z+1/2 
M i.41(x) -f M,(x)dx (6) 
z—1/2 
and 
z+1/2 
Ne41(x) -{ Ni(x)dx. (7) 
: z—1/2 
In view of 
atl/2 2sinu/2 | 
f e'“2dx = —______._ gius 
z—1/2 Uu 


we obtain (6) by an integration under the integral sign as follows 


z+1/2 1 0/2 sin u/2 k z+1/2 
[ M,(x)dx = - -f (= f etdx )du = My41(x). 
”) 2n —o Uu z—1/2 


¥ z-1/2 





To prove (7) we notice that 


k—1 
z+1/2 zx z Ny 
N,(x)dx = sf N,(x)dx = sf 6* ——_—_—— dx 
ie ; 2 : —2 (k 7 1)! 
k-1 k 
* X4. x 
= aot f eepromns ig on BP a ob Wace. 
-o (k— 1)! k! 


This concludes the proof of Theorem 3. 
3.13. Explicit polynomial expressions for M;(x). The formula 


M ,(x) = = ia (8) 


will readily show that y= M,(x) represents a spline curve of order k. Indeed, if k is 
even, then 
k-1 
(x + n)+ 


is a II, and therefore also their linear combination (8). If kis odd the same conclusion 
holds because 
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k—1 


(x-+n — 4)4 


is a IIx. 
It also follows from (8) that the explicit polynomial expressions for 


(k — 1)!M;(x), 


in successive unit intervals, are identical with the successive partial sums of the ex- 


pansion 


k\é-1 BE\ / k k-1 k\k-1 
bF¥yk-l = (« + =) - ( )(« +—- 1) +---+(- p*(« — =) ; 
2 1 2 2 


an expression which incidentally vanishes identically, being the kth order difference 
of a polynomial of degree k—1. We thus get 


i h 
0 if xs-— 
2 
k\e ie k ai k 
(:+>) ar Sail aiae ille 
k\Ft 7k k kl i 
(«+ ) -( ) («+ -1) if —-—+1s*s - +2 
2 1 2 2 
(kh —1)!Malx) = dere ccc cee renee eee n eee ene teen nee eenes (9) 
} k k-1 k k-1 
(+3) -G)G+2-") + 
k ; k A 


+ 
| 
aii 
“— 
| 
N| > 
4 
ed 
= 
= 
| 
| 
A 
IIA 
to 


+ 


| > 


| dtl = 0 if 


For future reference we work out explicitly the cases k =1, 2, 3, and 4. The expansions 


éx® = 1-1 
bx = (x + 1) — 2x + (x — 1) (10) 
6342 = (x + B)? — 3(x + 3)? + B(x — 9)? — (x — §)? 
643 = (x + 2)? — 4(x + 1)? + 6x — 4(x — 1)? + (x — 2)3 
now furnish the following expressions 
0 if x<—} 
M(x) =<1 if —}<x<}3 (11) 
0 if 4 < 2, 
to which must be added Mi(+1/2)=1/2 as required by (1) for k=1. 
| 0 if xs-1 
x+1 # -14<s 
piaiainal «+1 if O<x< an 
| 0 if 1s x, 
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(0 if xs-#% 
| (1/2)(a + 3/2)? if -#sx5-} 
M3(x) = 4 (1/2)(x + 3/2)? — (3/2)(x+1/2)? if -—}525} (13) 
(1/2)(— « + 3/2)? if 4 $Sxd3 
(0 if #S x, 
0 if x<-2 
(1/6)(* + 2)3 # —2425-1 
siti we | (1/6)( + 2) — (4/6)(x + 1)8 , ~tewes re 
(1/6)(— x + 2)* — (4/6)(— x + 1)’ if Osx 
(1/6)(— x + 2) if is<zs2 
(0 if 2S % 


In deriving these expressions the expansions (10) were used up to the point from where 
the evenness of the functions M;(x) allowed us to complete their definition for all x 
by symmetry. 


M, (x) 








ros 
~ 
' 
~ 
°o 
— 


—— 
° ae 
~ 
oe 
nD 
& 


4 & 











3.14. Interpolation formula with M;(x) as basic function. The interpolation formula 


ie} 


F(x) = > »Mi(x — ») (15) 


v=—oo 


will play an important role in our subsequent work. We mention it at this place be- 
cause it contributes to our investigation of k-order spline curves. 
To the basic function 


L(x) = M;(x) 
corresponds the characteristic function 
2 sin u/2\* 
g(u) = — Fe (16) 


as seen by comparing III (1) with II (22). The characteristic function of the formula 
(15) for integral x =n is 
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~! 
N 


ox(u) = M,(0) + 2M;,(1) cos u + 2M,(2) cos2u+---. (17) 


The {¢:(u)} represent an interesting sequence of cosine polynomials which we will 
investigate more closely later in this paper. Here we mention without proof that 


1 = o:(u) = d2o(u) > --- > d(u) > ---> O (0 < u < 2m) (18) 


while, of course, @.(0) =1. Hence (15) is a smoothing interpolation formula of progres- 
sively increasing strength as k increases. We assemble the various properties of (15) 
in the form of a theorem. 


THEOREM 4. 


F(x) = D> y,Mi(x — v) (19) 
is a polynomial smoothing interpolation formula of degree k—1, class C*-* and span 
2s=k (see sections 2.21 and 3.13). It is exact for the degree 1 and preserves the degree 
k—1. The smoothing power of (19) increases progressively for increasing values of k. 

The exactness of (19) for the degree 1 and the preservation of the degree k—1 
follow by Theorem 2 (B and C). Indeed, by (16), g(w) —1 has a double zero for u=0 
while g(u) has zeros of order k for u=27n (n #0). Since the preservation of the degree 
k—1 implies the identities I1(18) and (19), the following corollary results. 


COROLLARY. Any given polynomial P;-1(x) of degree at most k—1 may be represented 
in the form 


00 


P,. (x) = pe ¥nM (x —_ n) (20) 
where \yn} are the ordinates of some other suitably chosen polynomial of the same degree 
as P,1. This representation 1s unique. 

3.15. The analytic representation of spline curves of order k. We know that if } yn} 
is an arbitrary sequence of ordinates, then our interpolation formula 
F(x) = > V¥nM ,(x — n) (21) 
represents a spline curve of order k. This is true because all M(x —m) are such curves. 
The following question arises: Let F(x) be a given Il,; can we always represent it in the 
form (21) for an appropriate sequence \ yn} ? 
This question is answered affirmatively by the following theorem. 
THEOREM 5. Any spline curve I,(x) may be represented in one and only one way in 
the form 
II,(x) = > B ¥nM (x — n) (22) 
for appropriate values of the coefficients yn. There are no convergence difficulties since 
M;(x) vanishes for | x| >k/2. Thus (22) represents a I, for arbitrary \y,} and represents 
the most general one. 


In order to prove this theorem we return to the interpolation formula (21) and 
differentiate it repeatedly. By (8) we have 
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d 1 = ’ 1 = 
NiGes —_.__... 4, af ———- 4, 2 ae 
dx (k — 1)! (k — 2)! 
and repeating we get 
My (x) =8MiA(x) (OS v5 k—1). (23) 


From (21) and (23) we obtain by partial summation 


F'(x) = ps VrnOM (x — n) = 7. 6Yn41/2 Mia(x — n — 3) 


n n 
or 


F'(x + 4) = DO byn41/2Mu_s(x — n). (24) 


If k>2, this formal rule of differentiation of a spline curve may now again be applied 
to (24) with the result 
F"(x +1) = p 5° yn41M p_2(x% — n) 


a F(x) = ps 5°y,M;_2(x — n). 
Generally for OSvsSk-1 
b ™ 6’y,Mi.(x — n) if v is even, 
F(x) = \ (25) 


| >. 5’Vni1/2Mi(x —n— 43) if v is odd. 


This result may be stated as follows: The vth derivative of the spline curve (21) may be 
obtained directly by applying the same interpolation formula (21) with k—v, rather than 
k, to the sequence of the vth central differences 5’y properly centered according to the parity 
of v. In particularly: F*-®(x) is obtained by interpolating linearly among the 5*-*y. 
F%-)(x) is a step function whose successive values agree with those of the correspond- 
ing 6*—ly, 

Now let F(x) be a given II,. We are to show the existence of a sequence {yn} such 
that (21) holds identically. Suppose for the moment that such y, have been found 
which do make (21) hold. Then by (25) for y=k—1 we have 


be-ly, for n—i <x<n+ if kis odd. 
FURHD(g) = { 3 . (26) 


b*yauiye for n<x<cnt+il if k is even. 


In either case the successive constant values of the step function F“— (x) determine 
uniquely the values of the differences of order k—1 of the sequence of the as yet un- 
known coefficients y,. These differences in turn determine the coefficients y, uniquely 
up to an additive sequence of vanishing differences of order k—1. Let 5, be one se- 
quence such that 


5 Fnsiy2 = FO-V(n + 4) if Ris even 


or else 


b*15, = FO-D(n) if k is odd. 
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Consider the k-order spline curve 


F(x) 


I 


D> 5nMi(x — n) . 


and let 

R(x) = F(x) — F(x). 
From the way F(x) was defined it is clear that F“-)(x) and F@-»(x) agree in their 
successive unit intervals of constancy. Hence R(x) is a II, whose various polynomial 
arcs are of degree k—2 or lower. Therefore R(x) is identical with a polynomial of de- 
gree k—2. As such it allows of a representation of the form (21) in view of our Corol- 
lary of section 3.14. Therefore also 


F(x) = F(x) + R(x) 


may be represented by our formula (21). 
The unicity of the representation (21) is readily established. Indeed two different 


such representations would imply a representation of zero 


0= >> ynMi(x — n) 


without all y, vanishing. However the 6*~y all vanish, and our conclusion would con- 
tradict the uniqueness of the representation (20) of polynomials. 
A simple example might illustrate our proof of Theorem 5. Let us find the repre- 


sentation of the spline curve of order 4 


3 


1 
F(x) = Pn X4. 


By (26) we have 





s Fr" rs 1) ( 4 1)° { if n = 0 
a = n+%4)=(n+}3),= b 
aniads ‘ = 0 ft «<0. 
A sequence having these third differences is 
(O if n<0 
Vn =i f(n+1 n(n? — 1 
4 ( )- ( — tf a2 0. 
IX 3 6 
Hence 
ty, = Zz n(n? — 1)M,(x — n), (27) 
n=O 


with the possibility still open that both member might differ by a third degree poly- 
nomial. However this possibility is excluded by the remark that both sides vanish 
identically for x <0. Incidentally, (27) implies 


— (— x): = > n(n? — 1)M,4(x — n) 


n=—co 
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and by addition and subtraction of the two relations we get the identities 


= Do n(n?—1)Mx—n), — | x/*= DO | n(n? — 1)| Max — 0). 


In later applications we shall frequently operate with polygonal lines F(x) of de- 
gree k—1 not having continuity properties as strong as a IIx. Thus Jenkins’ L(x) 
defined by II(11) is of degree 4 and class C’’. Let F(x) be a polygonal line of degree 
k—1, having vertices at integral point x=, and being for all real x a function of 
class C* (—1Sy<k—2). We certainly obtain a curve of degree k—1 and class C* by 
addition of spline curves. 

F(x) = Wye + Mays +--+ + Tk, (28) 
where II, stands for II,(x) or II,(x+1/2), according to whether v is even or odd. 


THEOREM 6. Any given polygonal line F(x) of degree k—1 and class C* may be repre- 
sented as a sum (28) of k—y—1 appropriate spline curves of orders u+2, u+3,---,k. 


This theorem is a corollary to Theorem 5. Indeed, F“*+”(x) may have certain dis- 
continuities. We determine a II,;2 having the same discontinuities in its (u+1)st 
derivative. Then 

F(x) — Tl.+2 
is of degree k—1 and class C*+!. Proceeding in this way the theorem is readily estab- 
lished. 

Substituting for the II, in (28) their expressions in terms of the M, we obtain an 
explicit (parametric) representation of such polygonal lines. Thus Jenkins’ function 
II(11) may be represented as 


Liz) = 4M 4(x) + 3M a(x ot 1) of 4M,(x ~—- 1) —_ 2M;(x + 3) ~~ 2M;(x - 3). 


At a glance we recognize a curve of degree 4, class C’’, span s=6. 

3.16. A summation property of spline curves. The degree of a polynomial is de- 
creased or increased by one unit if we difference or else sum the polynomial. Not so 
with our spline curves II,. Indeed let 


k(x) = >> yaMi(x — n) (29) 


be a given spline curve. Then 
T(x + 1) = >> yaMi(x +1 — mn) = D> ynyiMi(x — n) 
and subtracting (29) we have 
AIT, (x) = I(x + 1) — T(x) = DO AyaMi(x — n). (30) 


Hence AIl;(x) will in general be also a IIx. Now let the spline curve II,*(x) be given 
and we wish to find II;(x) such that 


All, (x) = 11;(x). (31) 
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By using Theorem 5 the solution is immediately found. Indeed let 


IT ( 2) = 2, yn M 4 (x — n). 


Now (29) will be a solution of (31) provided Ay, =y,*, for all n. This relation defines 
the y, up to an additive constant which appears as an arbitrary additive constant in 
the solution II;,(x). It is thus seen that the operation of differencing or summing 
spline curves (29) reduces to performing the same operations on the sequence of the 
coefficients {y,}. 

3.17. An interpretation of M;(x) in probability theory. In concluding our discussion 
of polynomial spline curves we mention briefly the following interpretation of M;(x). 
As seen from its values as given by III(11) it is clear that Mi(x) may be interpreted 
as the probability density function of the error committed on a random real variable 
x, if that variable is rounded off to its nearest integral value. Now III(1) shows that 
the characteristic function of M;(x) is the kth power of the characteristic function of 
M(x). From a known proposition in probability theory we may conclude that M;(x) 
is the density distribution function of the error committed on the sum 


M+ tet: +s + % 
of k statistically independent real random variables x1, +--+ , xx, tf each variable is re- 
placed by its nearest integral value. 
This interpretation, otherwise entirely irrelevant for our purpose, does make a 
few of the properties of M;(x) intuitively obvious, such as 
>0O for |x| > k/2 
=0 for |x| < k/2, 


f M,(x)dx = 1. 


—oO 


In concluding we note the identity 


z+1/2 zy+1/2 rR-y+1/2 ea) 
f dn f dx2:- f f(x)dx, = f M,(u — x)f(u)du. 
zr 2 zi—1/2 zk-1—1/2 — 


—1/2 


3.2. Analytic spline curves of order k. The polynomial spline curves IIk(x) de- 
scribed in section 3.1 will be shown to be sufficient for the derivation of polynomial 
approximations to equidistant data enjoying various desirable properties. These poly- 
nomial approximations will have any a priori assigned number of continuous deriva- 
tives. However, in order to obtain analytic approximations we shall now proceed to 
derive from our spline curves II;,(x) an analoguous family of analytic functions. 

To achieve this end we shall smooth out our II;,(x) by means of one-dimensional 
heat flow. Consider an infinite homogeneous bar (the x-axis) in which the tempera- 
ture at the point x at the time ¢ is denoted by F(x, t). We assume the flow of heat to 
be governed by the equation 


fone petals (32) 


If F(x) = F(x, 0) is given, i.e., the temperature distribution at the time ‘=0 is known, 
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then F(x, t) is determined by the following integral'* 





F(x, 4) = — 
V4 


1 0 
= f e~(u—2)*t FB (y, 0) du. (33) 


This result is easy to verify: by partial differentiations we find that F(x, t) as defined 
by (33) indeed satisfies the differential equation (32) while familiar arguments origi- 
nated by Weierstrass will show that (33) implies 

lim F(x, 4) = F(x, 0), 

t++0 
provided F(x, 0) is continuous and, e.g., bounded. 

The solution of the problem of finding F(x, #), if F(x, 0) is given, is especially 

simple in the case when F(x, 0) is defined by a Fourier integral 


1 C) 
F(x, 0) = —f v(u)e*du. (34) 
2rd _«. 
Indeed, in this case we find 


1 «o 
F(z, = —{ e~ (ul) *Y(y)eiuzdu, (35) 
2x J _« 


Notice that the temperature ¢ enters only in the additional exponential factor. This 
can be proved in two ways, either by substituting (34) into (33) or else by verifying 
directly that (35) satisfies the differential equation (32). Obviously (35) reduces to 
(34) for t=0, as it should. 

We may (and wish to) think of F(x, ¢), for a fixed t>0, as a smoothed version of 
F(x) = F(x; 0). In fact F(x, t) is analytic and regular for all real or complex values of x 
if ¥(u) is, e.g., bounded. 

If we now apply this heat-flow transformation to our basic k-order spline curve 


1 */2 sin u/2\* | 
M(x) = f (- atti =) e™=du (36) 
2 J _« u 


we obtain by (35) its smoothed version 


1 “ 2/2 sin u/2\* | 
M(x, #) = — f e~ #(u/2) Gane edu. (37) 
2r J _« u 
Obviously 
M,(x, 0) = M;(x). (38) 


12 See H. S. Carslaw, Mathematical theory of the conduction of heat, Dover Publications, New York, 
1945, Chapter III, Section 16. Certain smoothing properties of heat flow were already noticed by Ch. 
Sturm in 1886. See in this connection G. Pélya, Qualitatives tiber Warmeausgleich, Z. angew. Math. u. 
Mech. 13, 125-128 (1933). It should be mentioned here that Weierstrass derived his famous approxima- 
tion theorem by means of the integral (33). Finally see E. Czuber, Wahrscheinlichkeitsrechnung, vol. I, 
Leipzig-Berlin, 1924, pp. 417-418, for a brief sketch of a method of using (33) to derive analytic approxi- 


mations to given data, 














78 I. J. SCHOENBERG [Vol. IV, No. 1 


The graph of y= M;(x, t) (t>0) is a bell-shaped curve which dampens out very fast. 
Later we shall learn how to compute its values very accurately. Here we mention that 


1 
0 < Mi(x, ) << —ae HD" Of og B R/2. (39) 
Vat 


Also, the recurrence relation III(6) generalizes so that Mi4:(x, ¢) is obtained by the 
averaging operation 


z+1/2 
M x41(x, t) = f M(x, t)dx. (40) 
z—1/2 
If now 
I(x) = >> yaMi(x — n) (41) 


is a spline curve of order k, then its heat-flow transform is 


io) 


II,(x, 4) = >> yaMi(x — 0, 2). (42) 
The graph of this function may be called an analytic spline curve of order k. We notice 
that the series (42) fails to converge only if y, increases very fast with | |. 
Summarizing we see that the curve (42) arises from the step function 


y y¥nM (x — n) 


n 


by k—1 successive applications of the averaging operation 


z+1/2 
f ( )dx 
z—1/2 


followed by the heat-flow transformation during a time interval t. The order of ap- 
plication of these k operations is of course irrelevant. 

The remaining parts of the paper are devoted to the problem of utilizing the two 
families of curves (41) and (42) for the purpose of approximating given equidistant 
data. 

IV. A FIRST CLASS OF ANALYTIC INTERPOLATION FORMULAE 


We shall now use the analytic spline curves III(42) to obtain 

1. A smoothing interpolation formula which is exact for the degree 1 only. 

2. An ordinary interpolation formula exact for the degree k—1. 

3. A smoothing interpolation formula depending on a positive parameter e, of 
strictly increasing smoothing power as ¢€ increases. For e=0 this formula reduces to 2, 
while for €= © it is identical with 1. 

4.1. A smoothing interpolation formula exact for the degree 1. A comparison of 
the formulae III(15) and III(42) immediately suggests the following analytic inter- 
polation formula 


x 


F(x) = D0 yaMi(x — n, 0) (1) 


n=—o 
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where M;(x, ¢) is defined by III(37). For the sequel we shall use the following notation 


,/2 sin u/2\* 
Valu, #) = e-tuin Ea) , (2) 
u 
‘2 sin u/2\* 
vi(u) = Vi(u, 0) - GS) ’ (3) 
u 
in terms of which III(37) becomes 
1 oo 
M(x, t) = —f ¥i(u, te™*du. (4) 
24 J _. 
The characteristic function of the smoothing formula (1) for integral x is 
ox(u, 1) = >> M,(n, 2) cos nu. (5) 
The general relation II(30) furnishes the following equivalent expression 
ox(u, t) = DY valu + 2a, 2). (6) 


The properties III(18) for the case =0 generalize for t>0 as follows 
1 > di(u, t) > g2o(u, 1) > --- > d(u, >--- >0, (0<u<2nr,t>0). (7) 


Moreover, for each fixed u, 0<u<2m, $:(u, t) is strictly decreasing as ¢ increases. 
These last two properties we state without proving them here. 

The arguments which lead to Theorem 4 now allow us to state the following theo- 
rem. 

THEOREM 7. For t>0, R=1,2,---, 

F(x) = }0 yaMi(x — 1, t) (8) 
is an analytic smoothing interpolation formula which is exact for the degree 1 and pre- 
serves the degree k—1. The smoothing power of (8) increases whenever either k or t is in- 
creased. 

4.2. An ordinary interpolation formula exact for the degree k—1. We shall now 
use the important property (7) to the effect that the periodic function ¢$;(u, ¢) is posi- 
tive for real u. This allows us to define the basic function 
1 = vi.(u, t) 


L,(x, tj) = — 
As, f 2x J _. o:(u, t) 


e“zdu (9) 





whose characteristic function is 


Vi(u, t) 
g(u) = Glad ° (10) 


Our Theorem 2 of Chapter II will readily yield the following result. 


THEOREM 8. For t20, k=1, 2,3,-:- 
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F(x) = > VnLi(x — n, b) (11) 


n=—oo 
is an ordinary interpolation formula which is exact for the degree k—1. 


Firstly we realize by (10) and (2) that the conditions II(33), (34) of Theorem 2 
are verified. Therefore (11) preserves the degree k—1. Secondly we have by (6), (10), 
and I1I(30) 


ital 
Se iA sealant 2 vi(u + 2nv, t) = 1. 
> 6b (a + 2a, 2) ox(u, t) “> 
Thus II(31) holds and (11) is therefore an ordinary interpolation formula. This con- 
cludes the proof of our theorem. Indeed, by the remark following Theorem 2 our 
formula (11) must also be exact for the degree k—1. 
We also mention without further details the following two limiting relations 


iat vi(u, 2) i vi(u, t) = { if | «| <T (12) 


i = lim ——— = ‘ 
t+o ,(u, t) hw z(u, t) Q if | u| > z. 


They show, in view of the integral representation (9), that our present basic function 
L(x, t) converges towards the basic function II (6) of the original cardinal series when- 
ever either & or else ¢ tends to infinity. 

4.3. A family of smoothing interpolation formulae depending on a smoothing 
parameter ¢. In section 4.1 we have derived the smoothing interpolation formula (8) 
in the derivation of which no attempt was made to compromise between smoothness 
of results and goodness of fit. Such a compromise is afforded by the following basic 
function 





o(u) = >> g(u + 2m) = 


IIA 


< wo) (13) 


1 * (u, t 
f ch ss Vi(u, the*du (0 


L,(x, t, 6) = — ketene 
' —o € + i(U, t)? 


2r 
which depends also on the smoothing parameter e. The corresponding interpolation 
formula 


eo 


F(x) = ps ynLi.(x — n, t, €) (14) 


n=—oo 


includes our previous formula (8) and (11) as special cases. Indeed by (13) ,(4) and 
(9) we find 


L(x, t, 0) = L,(x, t), (15) 
L(x, t, ©) = M;,(, 2). (16) 


Let us now investigate the characteristic function of the smoothing formula (14) 
for integral x. By (13) and II(30) this characteristic function is (¢;(u, t) is periodic!): 
e+ ¢:(u, t) 
o.(u, t, €) = ——————_ (u + 2/2, t). 
k ) ot ble, 1)? 2 vl / ) 


This and (6) give 
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blu, t, €-) = (epu(u, t) + ulus, t)*)/(e + be(u, t)*). (17) 
On the other hand we have by (7) the inequalities 
0 < g(u, 4) <1, (0 <u < 2n,t> 0). (18) 


Now (17), (18) imply 
0 < ¢gi(u, t, 6) < 1, (0<u<2r,t>0,€> 0), (19) 


and therefore (14) is a smoothing interpolation formula in the sense of section 2.21 b. 
Moreover, we see from (17) that for fixed ¢t and u (0<u<2r) ¢;(u, t, €) decreases 
monotonically from 
o(u,t,0) = 1 to gi(u,t, ©) = d(u, d) 

as € varies from e=0 toe= ©. Finally (14) is exact for the degree 1 and preserves the 
degree k—1 for the same reason as mentioned in the case of (8). 

We summarize these properties in the following theorem." 

THEOREM 9. For t20, R=1,2,---, 


2) 


F(x) = >> yLi(x — n, t, ©), (20) 


of basic function (13), is a smoothing interpolation formula which is exact for the degree 1 
and preserves the degree k—1. For €=0Q (20) is identical with the ordinary interpolation 
formula (11). For increasing values of € it increases in smoothing power until for €= « 
(20) is identical with our smoothing formula (8). 


4.31. A property of the derivatives of the approximation F(x). Let the given data 
\y¥n{ Satisfy the additional condition 


> | yn| < @. (21) 


n=—ow 


Assume also ¢>0, e€>0. We know that the sequence { F(n) } obtained by (20) is 
smoother than {| y,j in the sense of our discussion in I Section 1.1. However, it appears 


'8 The method used in deriving Theorems 8 and 9 obviously generalizes as follows. Let 
1 a2 
L(x) = — | g(uje' du 
2r Jw 
be a basic function. Let the periodic function 


o(u) = > g(u + 2xv) 


satisfy the inequality 
0 < ¢(u) £1. 


1 = g(u) 
loc — pte al turd 
L(x) 2x J. o(u) iin 


Then 


is the basic function of an ordinary interpolation formula. Moreover 
1 ¢* €+9(u) ' 
L(x, «) = — —— —— g(u)e*du 


Qed _~. e+ o(u)? 


gives rise to a family of smoothing formulae of increasing strength, as € increases. 
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to be of some interest to discuss here the smoothness of the function F(x), rather 
than that of the sequence { F(n) i, as a function of the smoothing parameter e. In 
this connection we prove the following 


THEOREM 10. Denote by F(x, €) the approximation (20) so as to indicate its depend- 
ence on €. The condition (21) insures the convergence of the integrals of the squares of 
the derivatives 


[em €))*dx, (m = 0,1, 2,---). (22) 


Also each of these integrals is a monotonically decreasing function of € in the range 
O<e< om, 
Indeed, let 


T(u) = >) ynei™ (23) 


be the characteristic function of the sequence {y,}. For convenience we define 


€ + ox(u, t) 








Q,(u, t, 6) = < sia 0 . (24) 
Then (13) becomes 
L.(x, t, 6) = —f Q,(u, t, eha(u, De~™*du. (25) 
By substitution of (25) into (20) we obtain 
F(x, «) = —f- Q,(u, t, ebWx(u, t) T(u)e~™*du. (26) 


We may evidently differentiate under the integral sign obtaining 


Se 
F(™)(z, 6) = — { Q,(u, t, vx(u, 2)T(u)(— tu) ™e—**du. 

2r J 0 
This formula exhibits the Fourier transform of F(x, €). We now use the analogue 
of the Parseval relation for Fourier integrals finally obtaining 


24? ™ du. (27) 





f em, e)*dx = — [out t, e)Wi(u t))?| T(u) 
= oe me 9% ’ 


This relation establishes our theorem. Indeed, by (24) the behaviour of the function 
0,(u, t, €) of period of 27, is as follows: Q;(0, t, €) =1, while for each fixed u, 0<u<2z, 
it decreases from 

Q,(u, t,0) = 1/dz(u, 4) to Q,(u, t, ©) = 1, 
as € increases from e€=0 toe= ~. 


The discussion of I section 1.2 concerning the smoothing of a finite table also has 
an analogue concerning the derivatives of the approximations F(x, e). We state the 
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result without further details. We assume the concrete situation of I section 1.2 where 
a finite table was extended to an infinite table by constant third differences at each 
end. To this extended table we apply our formula (20) with such a value of k which 
will insure that the formula (20) preserves cubics, i.e., 24. Then we can prove that 
the integrals (22) converge for m=4, 5, 6,- ++ and represent decreasing functions 
of «. 

4.32. Formula (20) as applied to subtabulation. Our formula (20) is excellently 
suited for the systematic interpolation, or subtabulation of given ordinates y,. It is 
less suited for interpolation. The reason is obvious: For subtabulation to tenths we 
need only a table of the basic function L(x, t, €) for the step h=0.1 only, while inter- 
polation would call for a much more elaborate table of this function. 

The following transformation of the formula (20), in terms of the function 
M(x, t), is of importance for numerical applications. First of all we expand the even 
periodic function (24) in a Fourier (cosine series) 
e+ ox(u, 2) = (k) 


p iu = (k) 
2,(u, t, -) = — ——— = w, (t,E&e = >, w, (t, 6) cos mm. (28) 
€E + bx ( ut, $)? v=—oo y=—w 





Substituting this expansion into (25) we get 


oo ’ 1 -) 
Li(z,tdo=- iw ¢o— f Vi(u, te-=)dy, 
Td 


rm—00 2 
which in view of (4) becomes 
Li(x,t,< = > wo, Ut, \Mi(x — », b). (29) 


This formula expresses the basic function L(x, t, €) in terms of the M;(x). If we sub- 
stitute this expansion into our formula (20) we obtain 


F(x) = Zz VnLi(% — n, t, €) - p » yao (ty -)Mi(x =e t) 


and replacing v by y—n 
F(x) = DY yawr-n(t, €-)Mi(x — », 0). 


A first summation by n introduces the sums 


fy = >> yuorn(ty €) (30) 


in terms of which our last expression becomes 
F(x) = >> f,Mi(x — », 2). (31) 


The pair of relations (30) and (31) is equivalent to (20) and represents its practical 
form. The reason for this is that the basic function M;(x, 4) dampens out like 
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exp(—x?) (see III (39)) while Li(x, t, €), e< ©, dampens out like exp(—x) only. We 
notice incidentally that (31) is identical with our formula (8), to be applied to the new 


computed ordinates {f,} given by (30). 
Frequently we require also tables of the derivatives F’(x) and F’’(x), of the ap- 
proximation F(x). These are then computed by the formulae 


F'(x) = >> f,Mi (x — », 0), (31’) 


> fMik' (x — », t), (31”) 


¥ 


F’'(x) 


from corresponding tables of Mf and My’. 

4.33. The least squares origin of formula (20). We want to sketch briefly the genesis 
of our formula (20). Let the sequence {yn} be given and consider the spline curve 
F(x) given by (31), where the coefficients {f.} are as yet unknown. If we try to deter- 
mine these unknowns by the requirement that F(x) should interpolate strictly the 
given ordinates yp, 1.e., 

F(n) = Yn (jo = 6 2 1, + 2,+*>), (32) 


we obtain an infinite system of linear equations in the unknown f,, the solution of 
which was found to be given by 


fe = Y 4y05-2(t, 0). 


n=—oo 


This leads to our ordinary interpolation formula (11). 

In view of Whittaker’s well known smoothing method it seemed natural to proceed 
now as follows: Let € be a given positive number. To determine the unknown coefficients f, 
of (31) as solutions of the following minimal problem, we set 


bt (F(n) — yn)? + €: y (fn — yp)? = minimum. (33) 


For e=0 the solution is of course identical with the solution of the ordinary inter- 
polation problem (32). For ¢€= the solution is obviously f,=¥y, in which case (31) 
reduces to (8). For0 Se @ a system of “normal equations” arises whose solution is 
found to be given by (30). The explicit solution of these normal equations (matrix 
inversion) is performed by the numerical determination of the cosine coefficients w, 
of the expansion (28) for each given set of values of k, ¢ and e. 


V. THE COMPUTATION OF THE TABLES I, I, Ill 


In this last chapter we shall discuss the methods used in the computation of our 
tables which allow us to use our formulae (30), (31) or 


fn = DX yoonv(t, 0), (1) 
F(x) = >> faMi(x — n, 0), (2) 


for subtabulation to tenths. 
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5.1. The computation of the function M;(x, t) and its derivatives. This function 
is defined (see III(37)) by the integral 


. 2 , (2 sin u/2\* | 
M(x, t) = —f e~ t(ul2)" {| —___ }_ giuzdy, (t> 0). (3) 
209 u 


For the special case of t=0 and k=1, 2,3, - - - we found previously the explicit poly- 
nomial expressions III(3). It now so happens that for our present case of ¢>0 (3) al- 
lows us to define our function also for k=0 as 


1 0 
M(x, t) = —f e~ *(ul2) eiuzdy, 
2a J 
This last integral is known to be identical to 
1 2 
M(x, t) Ss ” iad it, (4) 
V rt 


The recurrence relation III(40) shows that (3) is obtained from (4) by repeating k 
times the averaging operation 


z+1/2 z 
f or sf ’ 
z—1/2 —x 
The result, however, is not changed if we perform all k integrations first to be fol- 


lowed by the operation 6* of kth central differencing. This proves the following result: 
If we define a sequence of functions g;(x, t) by 





1 " 
£o(x, t) aaa ~— re (5) 
V tt 
and the recurrence relation 
gu(x, t) -{ gia(x, idx (k= 1,2,3,---), (6) 
then 
M(x, t) = 6*g,(x, 2). (7) 


This relation reduces the problem to the problem of computing the repeated in- 
tegral g,(x, t) of the error function (5). This we do as follows. It is easy to prove by 
induction or otherwise that (5) and (6) imply 


1 1 z 4 
g.(2, 2) = ———_ — f{ (x — u)*'e-“ ! du. (8) 
(kR-—1)! VariJd_. 
With x—u=v this becomes 
1 1 e 
g(x t) is: ates — | e7 (2-9)? tyk—-Idy, (9) 
pic (k—1)! VartJo 


By differentiating this with respect to x we get 
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1 1 ” 


‘(z,) = ——_ -—— 
ge (%, 1) (k—1)! Vxtdo 


eW (2-H tyk-1(— 2xt-l 4 Q0t-) do 
2x 2k 
=— = g(x, t) + .y £41(%, 2) 
and therefore the recurrence relation 


t x 
geval, ) = ai (x, 0) + = gu(a _ 00ee.<s3 (10) 


which allows us to compute the successive values g,(x, t) by the operation of differen- 
tiation rather than integration. Indeed from (5) and (6) for k=1 we get 


1 z 
gi(x, t) = ; f e~*'/tdx, (11) 


V rt 





Now (11) and (10) for k=1 will give 








t 1 ‘. 1 *.. 
go(x, 4) = — ez lt 4 x f e~* |'dx 
2 Vat V tid —« 


from which g;(x, ¢) is readily determined. 
This progressive computation is greatly simplified if we realize that gi(x, ¢) will 
appear as an expression of the form 


Bex, t) = Pr(x, t)go(x, t) + Q(x, dgrl(x, 2), (12) 


where P;, Q; are polynomials in x and ¢, while go(x, ¢) is the error function (5) and 
gi(x, ¢) is the error integral (11). Substituting (12) into (10) we find 


t t x 
Beri(x, ) = ok (Pk + Ox) go(x, t) + (x + x 0) 81 (x, t). 
On comparing with (12) for k+1, rather than k, we obtain the recurrence relations 


t 
Pix = > (Px + Q:) 


Qui = Of += 0 (k =1,2,---). (13) 
Since P;=0, Qi:i=1 we readily obtain the following explicit expressions 
P, = t/2, Q2 = x, 
P; = tx/4, QO; = t/4 + x*/2, 
P, = t(t+ x*)/12, QO, = x(3¢ + 2x?)/12, (14) 
Ps; = tx(St + 2x?)/96, Os = (37 + 12tx? + 4x4) /96, 


x(15é2 + 20tx? + 4x4) /480. 


P, = t(4t? + 9tx? + 2x4) /480, Os 


Excellent tables of the probability function (5) and its integral (11) are now available. 
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By means of these tables the formulae (12) and (14) allow us to compute readily the 
function g;(x, ¢).14 It seems worth while to point out that the relation (7) goes over 
into III(8) if +0. Indeed by an obvious change of variable we see that (11) becomes 
1 z/ vt a 
£i(x, t) ne = e~“ du 
V TH -x 
and therefore 


0 
lim gi(x, t) = x4. 
t+ +0 


Now by induction we prove by (6), on letting +0, that 


1 k-1 


lim g(x, t) _ + 


aoe: 
t++0 (k = 1)! 


which proves our last statement by continuity. 
The computation of the derivatives of M;(x, t) is immediately settled by the rela- 


tion 


My (x, t) = 6 gu-0(x, 2), (15) 


which is implied by (6). 

5.2. The computation of the cosine coefficients w\") (t, «). By IV(5) and IV(28) 
we can see that the problem consists in computing the values of the coefficients of 
the cosine expansion of the function 

Q;(u, t, 6) = stone => on (b, €) COS nu, (16) 
e+ di(u, 2)? = nm —eo 


where the even periodic function ¢;(u, 4) is defined by its cosine expansion 


o:(u, t) = > M ,(n, t) cos nu. (17) 


n=—o 


14 | owe to D. H. Lehmer the reference to the functions Hh,(x) defined by 
Hho(x) = f e-**/2dx, Hhp(x) = f Hhn_-1(x)dx. 


Tables of these functions were published by J. R. Airey as Tables XV, Group IV, of the Mathematical 
Tables of the British Association for the Advancement of Science. The relation between our gx(x, ¢) and 
these new functions is 


1 
ge(x, t) = —=— (t/2)* 2A a(— xV2/t). 
V2 


This relation, for k=4, t=}, would readily allow us to compute our Table I by means of Airey’s tables 
of Hh;, Hhs and Hh. However, for other sets of values of k and ¢, such as k=8,t=}, which are needed 
for other purposes, the range of x in Airey’s table becomes insufficient. In this case tables of M,(x, t) and 
its derivatives are computed by our formula (12) and the excellent Tables of Probability Functions, vol. I 
(1941), vol. II (1942), prepared by the Mathematical Tables Project under the direction of A. N. Lowan. 
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The coefficients /;(n, t) of this extremely fast convergent series are readily computed 
to 10 decimal places. The obvious procedure would be to compute from (17) a table 
of ¢;(u, ¢), then compute a similar table of Q;(u, ¢, €) which is then to be used in 
computing w, by some method of numerical harmonic analysis. It would be hard to 
achieve accuracy by this method and for this reason we proceeded differently. It should 
be born in mind that the cosine expansion of the denominator of (16) is readily ob- 
tained by the simple operation of multiplication of Fourier series. The only trouble- 
some part is the computation of the expansion 


1 i) 
—_——— = > &, cos nu, (18) 
et gi(u,? —. 





i.e., the reciprocation of a given cosine series. This was done as follows. The above- 
mentioned method of a 24-ordinate harmonic analysis scheme was used for obtaining 
values of the c,’s accurate to 4-5 decimal places. These values were then improved 
to values accurate to 8 or 9 places by an iteration method developed by H. A. Rade- 
macher and the author. This method is closely related to the method recommended 
by H. Hotelling for the reciprocation of ordinary matrices and will be described else- 
where. 

In concluding this paper we want to point out two special cases of our ordinary 
interpolation formula (11), or (20) for €=0, which are of mathematical interest. We 
mention first the case of k=0, ¢>0. This corresponds to interpolating our ordinates 
yn by means of a function F(x) as described by the formula (8) of the Introduction. 
Although, as remarked there, the resulting interpolation formula is useless for practi- 
cal purposes, it has the remarkable feature that the expansion coefficients w(t, 0) 
of (16) may be obtained explicitly. Indeed the function ¢o(u, ¢) reduces to a Theta 
function which is a regular and uniform function of 


= piu 
Z=€ 


with singularities only at z=0 and z= «. The simple zeros of this function are real, 
negative and form a geometric progression. As a result we are able to find explicitly 
the decomposition in partial fractions of the reciprocal 


1 /do(u, t) ° 


The expansion of these partial fractions into geometric power series furnishes ex- 
plicitly the Laurent expansion in powers of z and therefore also the cosine expansion 
(16). 
The second special case of interest is k>0, t=0. In this case our formula (11) re- 
duces to an ordinary polynomial interpolation formula of degree k—1 and class k—2. 
This does not contradict Mr. Greville’s statement (loc. cit. page 212) to the effect 
that such formulae do not exist. Indeed Mr. Greville considers only basic polynomial 
functions L(x) of finite span s only, while our basic L;(x, 0) are of infinite span. This 
case, which is of considerable interest, requires a more detailed investigation of the 
cosine polynomials ¢;(u, 0). We postpone this discussion to the second part B of 


this paper, 
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APPENDIX 


Description of the tables and their use for the analytic approximation 
of equidistant data. 


Tables I and II. In Table I we find the 8-place values of the even function 
17° ,./2sin #/2\* 
M(x) = —{ e* n(—*) cos uxdu (1) 
2x J _. u 


and its derivatives M’(x), M’’(x) for the step of Ax =0.1. The graph of M(x) is a bell- 
shaped curve and M(x) vanishes to 8 places for x S$ —4.3 and x24.3. We now define 
a function of period 27 by the cosine series 





¢(u) = M(0) + 2M(1) cos u + 2M(2) cos 2u+-:- (2) 
and expand in cosine series the following functions 
e+ $(u) 
———— = wo(e) + 2wi(e) cos u + 2uwe(e) cos2u+---, (3) 
e+ o(u)? 
where € is a non-negative parameter. Our Table II gives the 8-plane values of these 
coefficients for €=0, 0.1, 0.2, ---, 1.0. 


These tables may be used as follows to obtain an analytic approximation F(x) 
to our ordinates y,. We discuss first the case when F(x) is to interpolate the ordinates, 


in the usual sense, i.e., 

F(n) = Yn. (4) 
For this end we compute first from the sequence {y,} a new sequence of coefficients 
{fn} by means of the formula 
fr = +++ + yn—202(0) + yn—101(0) + ynwo(0) + Yn4101(0) + yn+2m2(0) + -°- (S) 


or 


te = z Yrn—(0), (5’) 


where Wp» =@_m. The analytic approximation of the ordinates y, is then given by 


F(x) = D faM(x — 1). (6) 


The values of F(x), x to even tenths, are now readily computed. Thus 


F(2.3) = f-1M(3.3) + foM (2.3) + fiM(1.3) + foM (0.3) 
+ faM(— 1+.3) + faM(— 2 +.3) + foM(— 3 + .3) + foM(— 4+ .3). 


The tabular values of M(x) are so arranged that all 8 values needed in this computa- 
tion are found in the fourth column headed x+.3. Generally, if the values of f, are 
written in a vertical column, we compute the values of F(m+v-10-') (v=0,1, -- -, 9) 
by matching the column of values of f, with the vth column of the table of M(x) in 
such a way that f» corresponds to the row for x =0. The products 


{,M(m — n + v-10—') 
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are then accumulated in the products counter of a desk computing machine. Also the 
values f, are best computed by (5’) in a similar way if the column of values of w,(0) 
is extended upwards by symmetry for negative values of n. 

From the tables of M’(x) and M’’(x) we may likewise compute tables of the de- 
rivatives of F(x) by 


F(x) = > feM r(x _ n). (7) 


A check of the computation of the coefficients f, is afforded by (4). Indeed the values 
F(n) computed by (6) should agree with the y, to about eight significant figures. 
The formula (6) is exact for cubics, i.e., if the y, are the ordinates of a polynomial 
of degree at most 3, then F(x) is identical with that polynomial. 
If the conditions (4) of strict interpolation are not required, then we have the 
possibility of obtaining an approximation F(x) which is such that the sequence 
| F(n)} is smoother than the given {yn}. The approximation F(x) is then given 


by the pair of formulae 


> 
| 


Dd yrn+(6), (8) 
> fa (x —n), (9) 


Il 


F(x) 


which are applied as above. The choice of the value of the smoothing parameter € 
depends on the amount of smoothing desired. The strongest smoothing afforded by 
our table is obtained for e=+. Then (3) shows that wo(o)=1, wi(o) =we(o) 
= --- =0. Thus (8) becomes f, =y, and (9) reduces to 


F(x) = >> yaM(x — n). (10) 


This formula is especially simple to apply. It should be remarked however that, if 
€>0, our formula (9) is exact only for linear functions and the same is true of (10). 

Table III. We may eliminate the coefficients f, between (8) and (9). In terms of 
the new even function 


L(x, 6) = >> wale)M(x — n), (11) 


our formulae (8), (9), then reduce to 


F(x) = B YnL(x — n, €). (12) 


Table III gives the values of L(x, €) and L’’(x, €) for e=0, 0.1, - - - , 1.0 for the step 
Ax =6.5. These may be used for subtabulation to halves in preference to (5), (6) or 
(8), (9). For subtabulation to fifths or tens, the use of formulae (8), (9) is preferable 
because of the slower damping of the function L(x, €). Even so, formula (12) and 
Table III allow us to estimate quickly how well F(x) approximates the y,. By (12) we 


have 
F(x) = >> yal’ (x — n, ©). (13) 


The table of L’’(x, €) then allows us to compute quickly a table of F’’(x) for the step 
Ax =0.5 or else only isolated values if such are needed. 
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Example of subtabulation to tenths. We consider the following fairly smooth se- 
quence of 64 ordinates y,: 














n Vn n yn || m | ym |] m | Yn || 1m | Yn ! n | Yn 
| 1] 1} | 
1 24614 12 | 25370 | 23 | 29290 | 34| 73820 || 45 | 82450 | 56 | 79962 
2 24644 13 | 25504 } 24 | 30160 || 35 | 77830 || 46| 82290 || 57 | 79698 
3 | 24680 14 | 25660 25 | 31320 || 36| 80240 || 47 82110 | 58 | 79431 
4 | 24723 is | 25850 || 26 | 32840 |) 37| 81660 || 48| si911 | 59] 79161 
5 | 24772 16 | 26080 || 27 | 34790 || 38] 82330 || 49 | 81699 || 60 | 78889 
6 | 24828 17 | 26350 || 28 | 37260 || 39 | 82680 || 50 a 61 | 78614 
7 24892 18 | 26660 || 29 | 40440 || 40 82840 || 51 | 81234 || 62 | 78338 
8 | 24966 19 | 27040 || 30 | 44750 | 41] 82830 || 52 | 80987 | 63 | 78060 
9 | 25048 20 | 27490 || 31 | 51120 || 42) 82780 || 53 | 80736 || 64| 77780 
10 | 25143 21 | 28010 || 32 | 59390 || 43 | 82700 | 54 | 80481 | 
| 67550 || 44] 82590 || 55 | 80223 | 














11 25250 22 


28600 


33 








The differences of the section of this table with which we will be concerned are as 











follows: 
n Yn A | At | a3 | At | as 
27 34790 | | | 
28 37260 2470 | 
29 | 40440 3180 | 710 | 
30 | 44750 | 4310 ‘| 1130 | 420 
31 | 51120 6370 2060 | 930 510 
32 59390 8270 | 1900 | — 160 —1090 —1600 
33 67550 | 8160 | -— 110 | 2010 | -1850 | — 760 
34 | 73820 | 6270 —1890 | —1780 | 230 | 2080 
35 | 77830 | 4010 | -2260 | — 370 1410 | 1180 
36 | 30240 2410 | —1600 | 660 1030 | = — 380 
37. | 81660 1420 - 990 | 60 | — so | 1080 
38 | 82330 670 — 750 200 | %—370 | -— 320 








We illustrate the case of strict interpolation, i.e., we use our Tables II for e=0. From 
our formula (5) and the values of w, as given in the column of Table II, with the head- 


ing €=0, we obtain the following coefficients. 





n fn 

27 34662 .222 
28 37031 .355 
29 40215 .195 
30 44060. 182 
31 50349 .304 
32 59490 .524 
33 68212 .510 
34 74566 .216 
35 78283 .074 
36 80460 .234 
37 81953 .811 


38 82356.888 
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From these values and our Table I of M(x) and M’’(x), we obtain by the formulae 
(6) and (7) the following tables of F(x) and F’’(x) with their differences. 


Table of the function F(x) and of its second derivative sles 






































x F(x) 4 | A? At F(x) | A | 2 a? | at 
31.0 | $1120.00 | 2117.97 | | | | 
31.1 | 51884.17 | 76417 | | | || 1966.48 | —15149 | | 
31.2 | $2667.97 | 78380 | 1963 | | 1787.44 | —17904 | —2755 | 
31.3 | 53469.63 | 80166 | 1786 | —177 | 1583.71 | —20373 | —2469 | 286 | 
31.4 | 54287.11 | 81748 | 1582 | —204 | —27 || 1359.15 | 22456 | —2083 | 386 | 100 
31.5 | 55118.17 | 83106 | 1358 | —224 | —20 | 1118.30 | —24085 | —1629 | 454 | 68 
31.6 | 55960.40 | 84223 | 1117 | —241 | —-17 || 866.08 | —25222 | —1137 | 492] 38 
31.7 | $6811.29 | 85089 | 866 | —251 | —10 | 607.04 | —25868 | — 646 | 491 | — 1 
31.8 | 57668.25 | 85696 | 607 | —259| — 8 || 346.89 | —26051 | ~< $63 | 463 | —28 
31.9 | 58528.68 | 86043 347 | —260 | — 1 || 88.63 | —25826 225 | 408 | —S55 
32.0 | 59390.00 | 86132 89 | —258 2 || — 163.98 | —25261 565 | 340] —68 
32.1 | 60249.69 | 85969 | — 163 | —252 6 || — 408.22 | —24424 837 | 272 | —68 
32.2 | 61105.30 | 85561 | — 408 | —245 | 7 {|| — 642.14 | —23392] 1032] 195 | —77 
32.3 | 61954.51 | 84921 | — 640 | —232| 13 || — g64.26| —22212| 1190] 148 | —47 
32.4 | 62795.08 | 84057 | — 864 | —224 8 || —1073.51 | —20925 1287 | 107 | —41 
32.5 | 63624.93 | 82985 | —1072 | —208] 16 || —1269.11 | —19560| 1365 | 78 | —29 
32.6 | 64442.10 | 81717 | —1268 | —196 | 12 || —1450.39 | —18128 | 1432] 67 | —11 
32.7 | 65244.77 | 80267 | —1450 | —182 | 14 } —1616.76 | —16637 | 1491] 59| - 8 
32.8 | 66031.30 | 78653 | —1614 | —164 18 } —1767.70 | —15094 | 1543| 52| —7 
32.9 | 66800.16 | 76886 | —1767 | —153 | 11 || —1902.77 | —13507 | 1587| 44| — 8 
33.0 | 67550.00 | 74984 | —1902 | —135 | 18 || —2021.68 | —11891 1616 | 29 | —15 
33.1 | 68279.64 | 72964 | —2020| —118| 17 | —2124.30 | —10262 | 1629] 13 | —16 
33.2 | 68988.05 | 70841 | —2123 | -103 | 15 | —2210.71 | — 8641 1621 | — 8 | —21 
33.3 | 69674.37 | 68632 | —2209| — 86/ 17 || —2281.13 | — 7042 | 1599 | —22 | —14 
33.4 | 70337.91 | 65354 | —2278| — 69 | 17 || —2335.91| — 5478| 1564] —35 | —13 
33.5 | 70978.07 | 64016 | —2338 | — 60 9 || —2375.46 | — 3955 1523 | —41 | — 6 
33.6 | 71594.50 | 61643 | —2373 | — 35 25 || —2400.17 | — 2471 1484 | -39| 2 
33.7 | 72186.94 | 59244 | —2399 | — 26 9 || —2410.41 | — 1024] 1447 | —37 | 2 
33.8 | 72755.29 | 56835 | —2409 | —10| 16 || —2406.55 386 | 1410 | —37| 0 
33.9 | 73299.58 | 54429 | —2406 | 3 13 || —2389.01 | 1754] 1368 | —42| — o 
34.0 | 73820.00 | 52042 | —2387 | 19] 16 || 2358.32] 3069} 1315 | —53 | — 











An inspection of these tables shows that they are very smooth and that they define 
F(x) and F’(x) to 7 significant figures by 4-point central interpolation. We have 
chosen on purpose an example for which it would be hard to obtain similar results by 
standard methods, if we are to maintain the forced accuracy requirement, and the 
same high degree of consistency between the function F(x) and its second derivative 
F’' (x). For purposes of comparison we show also the interpolated values F,(x) for the 
range x =31.6—32.5 obtained by the 10-point central interpolation method. 

On comparing with our table of F(x) we notice that 

F(x) < F(x) 


throughout this range, with the exception of the point x =32.0 where, of course, both 
values agree. The curve F,(x) has a corner at x =32. This is the typical discontinuity 
in the first derivative due to central interpolation methods (see the first paragraph 


of our Introduction). 
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x F.(x) | A A? A’ Af 
31.6 55959 .90 
31.7 56810.60 | 85070 | 
31.8 57667 .55 85695 | 625 | 
31.9 58528.22 | 86067 | 372 —253 
32.0 59390 .00 | 86178 } 111 —261 — 8 
32.1 60248 .72 85872 | — 306 —417 —156 
32.2 61103.61 | 85489 | — 383 | — 77 340 
32.3 61952.37 | 84876 — 613 —230 | —153 
32.4 62792.77 | 84040 — 836 —223 | 7 
32.5 


63622 .68 82991 — 1049 —213 10 


Notice that we needed 12 coefficients f, for the subtabulation of three panels. 
Each additional coefficient f, (n=39, 40, - - - ) allows the subtabulation of an addi- 
tional panel. 

It should be remarked that 53 ordinates y, enter into the computation of each 
coefficient f,. This is due to the slow rate of damping of the w,(€) for «=0. Thus for 
€=.1 (very moderate smoothing) only 35 ordinates y, are needed, for €=1.0 only 23, 
for €= © only 1. Concerning the important matter of dealing with the ends of a table 
see section 1.2 and the last paragraph of section 4.31. 


TABLE I: M,(x, t), Mi (x, }), Mi’ (x, #) for k=4, t=0.5, Ax=0.1. 
M(x, 1/2) 


x+.4 











x x+0 x+.1 | x+.2 x+.3 
4 00000004 .00000002 .00000001 | 
3 00011325 00005910 | .00002991 .00001467 .00000697 
2 01616917 .01105340 .00737858 | .00480621 | 00305258 
1 22597004 18940616 15590118 12596479 | .09986387 
0 51549499 51132566 | .49901141 47911917 | .45254731 
1 22597004 26483185 30499058 34523755 | 33420963 
-2 01616917 02311310 | 03230776 04418973 | .05917998 
3 00011325 00021062 | 00038032 00066726 | 00113822 
4 .00000004 .00000010 .00000026 .00000062 | .00000143 
x x+.5 x+.6 x+.7 x+.8 | x+.9 
3 .00000321 .00000143 .00000062 .00000026 | .00000010 
2 00188907 .00113822 | .000667 26 00038032 00021062 
1 07764689 .05917998 | .04418973 03230776 | .02311310 
0 .42046084 .38420963 .34523755 | .30499058 | 26483185 
1 .42046084 .45254731 AT9II917 | .49901141 | .51132566 
3 07764689 09986387 12596479 .15590118 | .18940616 
— 00188907 .00305258 | .00480621 .00737858 01105340 
4 00000321 .00000697 | 00001467 00002991 .00005910 
—5 00000001 .00000002 
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Mi (x, 1/2) 





























e | x+0 x+.1 | x+.2 x+.3 | x+.4 
| | 

4 — .00000039 | —.00000015 | —.00000006 — .00000002 | —.00000001 
3 = .00071955 | — .00039340 — .00020833 — .00010680 | —.00005298 
2 — .05961795 | —.04334506 i= .03071644 — .02120376 | —.01424920 
1 — .37860391 — .35140346 | —.31784825 — .28043885 = .24150489 
0 | .00000000 | —.08306134 | — .16227165 — .23406492 | — .29541674 
—1 .37860391 -39695855 -40419276 .39846265 | .37855467 
—2 .05961795 | .07996844 . 10465732 - 13368990 . 16673619 
-—3 | .00071955 .00127546 .00219236 .00365652 .00592117 
—4 .00000039 .00000096 | .00000229 .00000529 .00001179 

x x+.5 x+.6 | x+.7 | x+.8 | x+.9 
3 — .00002542 | —.00001179 ze -00000529 | .00000229 | — .00000096 
2 — .00931577 |} —,.00592117 — .00365652 | — .00219236 | —.00217546 
1 ee . 20306520 | —.16673619 | —.13368990 — .10465732 | —.07996844 
0 — .34404758 — .37855467 — .39846265 | — .40419276 | —.39695855 
—1 .34404758 .29541674 . 23406492 | . 16227165 -08306134 
—2 . 20306520 . 24150489 . 28043885 .31784825 -35140346 
~“s | .00931577 .01424920 .02120376 .03071644 .04334506 
—4 .00002542 | .00005298 } .00010680 .00020833 | .00039340 
| .00000006 .00000015 


=—§ | .00000001 .00000002 





Mit! (x, 1/2) 











x | x+.0 x+.1 x+.2 | x+.3 x+.4 





























| 
4 | 00000357 00000145 | — .00000056 .00000021 | — .00000008 
3 | — .00423106 .00243772 .00135797 .00073109 .00038024 
2 18251117 14368197 10978191 .08140988 .05858190 
.23181861 .30800376 | — .35890239 .38537940 .38991971 
0 | —.83712882 — .81763132 | —.76058819 — .67020231 — .55301267 
-1 23181861 13144694 .01013023 | — .12678241 — .27209706 
—2 18251117 .22494722 26885281 .31126005 34845442 
—3 |  .00423106 00710375 01154277 .01816077 .02768079 
-4 .00000357 .0000085 1 .00001955 .00004332 .00009259 
—5 | | 00000001 

| | 
x +5 | r+.6 | x+.7 | x+.8 x+.9 

| | 

4 | 00000003 | 00000001 ‘| 
3 | 00019097 | — .00009259 | —.00004332 | © .00001955 .00000851 
2 .04089359 | 02768079 | .01816077 01154277 .00710375 
1 .37617315 |  .34845442 | —.31126005 . 26885281 .22494722 
o | —.41725773 | —.27209706 | —.12678241 .01013023 13144694 
—1 | —.41725773 | —.55301267 | —.67020231 — .76058819 — .81763132 
—2 | 37617315 |  .38991971 .38537940 .35890239 .30800376 
-3 04089359 | 05858190 | .08140988 | 10978191 14368197 
—4 | 00019097 | 00038024 | 00073109 { .00135797 .00243772 
—5 .00000003 | 00000008 | .00000021 |  —.00000056 .00000145 
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TABLE II: w(t, €) for k=4, #=0.5, «=O (0.1) 1.0. 
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«=.0 | e=.1 e=.2 ex 3 e=.4 
3.50637741 | 1.61378653 1.39953009 1.30308904 1.24631521 
—1.84900618 | —.26890929 — .14132793 — .09293505 — .06784110 
87238793 | —.08981772 — 09332063 — .08242675 — .07223261 
— 40443570 | .07027891 04023694 .02480160 01624479 
.18693997 | —.02078617 | —.00397484 .00050538 .00184671 
— 08636451 | 00133949 | —.00223908 — .00188468 — .00132999 
03989615 | 00169234 .00099749 00037463 00010684 
— 01842978 | —.00088114 — .00010447_ | .00005298 00006541 
00851350 | -00019734 | —.00005372 — 00003828 — .00001761 
— 00393275 | 00001073 | 00002469 00000455 — .00000116 
00181670 | —.00002625 — .00000273 .00000174 .00000130 
— 00083921 | 00001022 — 00000129 — .00000070 — .00000014 
00038767 | —.00000156 00000061 .00000003 — .00000006 
— 00017908 | —.00000044 — .00000007 .00000004 .00000002 
00008272 | .00000036 — .00000003 | —.00000001 
— 00003821 | —.00000011 | 00000002 | 
.00001765 | .00000001 | 
— 00000815 00000001 | 
.00000377 
— .00000174 
.00000080 
— .00000037 | 
.00000017 | 
— .00000008 | 
.00000004 si 
— .00000002 
.00000001 
e=.5 e=.6 | e=.7 e=.8 | e=.9 | «=1.0 
20834767 18095463 | 1.16016154 1.14379093 | 1. 13054096 1.11958158 
05268720 | —.04264921 | —.03556878 | —.03034057 | —.02634263 | —.02319971 
.06387537 | —.05710538 | —.05156939 | —.04698035 | —.04312407 | — .03984269 
01107108 | .00773445 | 00547641 .00389078 | 00274451 .00189634 
00219079 00218246 | .00204868 | .00187682 00170183 00153738 
00091344 | —.00062674 | —.00043120 | —.00029659 | — .00020265 | —.00013620 
.00000351 | —.00004704 | — .00006163 | — .00006358 | —.00006018 | — .00005476 
.00005180 | .00003700 00002545 | 00001716 | .00001138 .00000739 
00000656 | —.00000140 |  .00000085 | = .00000171 | = .00000193 .00000187 
- 00000206 | —.00000175 | —.00000127 | —.00000087 | —.00000057 | — .00000036 
.00000064 | .00000025 | .00000006 | —.00000002 | —.00000005 | — .00000006 
00000003 | — .00000006 | ~— .00000006 | = .00000004 |  .00000002 |  — .00000001 
.00000001 | 


.00000004 = 


.00000002 : 
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TABLE III: L(x, t, €), Le’ (x, t, &) for k=4, #=0.5, e=0 (0.1) 1.0, Ax=0.5. 
Li(x, 1/2, © 

x e=.0 e=,1 | «= e=.3 e=4 
0.0 | — 1.00000000 70747935 | 65457028 .62707488 60947694 
0.5 62191163 53757743 | .51070485 .49509729 48452022 
1.0 .00000000 .20252568 | .22066478 .22681461 .22949350 
1.5 — 17291085 — 02061576 | .01285810 .02919893 03901344 
2.0 .00000000 —.06545791 | —.04840114 .03681939 — .02872086 
2.5 .07415615 — .02765903 | —.03096280 .02894823 — .02631330 
3.0 .00000000 00709183 | —.00340667 00711221 — .00850828 
3.5 — .03382251 01344008 | .00756627 .00392340 .00177090 
4.0 .00000000 .00401304 .00502862 .00410188 .00314843 
4.5 .01558996 — .00276042 .00041208 00121901 .00135036 
5.0 .00000000 — 00251219 — .00118870 .00038015 .00001138 
5.5 — .00719897 — 00027479 | —.00076318 .00054505 — .00033539 
6.0 .00000000 00065102 | —.00007596 .00021043 — .00019927 
6.5 | .00332530 .00042139 | .00019012 .00003152 — .00002867 
7.0 | .00000000 — .00000773 .00012316 .00007297 .00003257 
7.5 | —.00153608 — .00015286 .00000861 .00003207 .00002580 
8.0 | 00000000 | —.00006787 — .00002989 .00000086 .00000704 
8.5 | .00070958 | .00002025 | —.00001865 .00000923 — .00000251 
9.0 | 00000000 | .00003031 | —.00000162 .00000502 — .00000321 
9.5 | —.00032779 | .00000793 | 00000477 | .00000028 — .000001 20 

10.0 | .00000000 | — 00000573 .00000301 | .00000115 .00000010 

10.5 | 00015142 | — .00000566 .00000017 | .0000007 2 .00000036 

11.0 .00000000 — .00000083 —.00000075 | .0000001 1 .00000018 

11.5 | —.00006995 | .00000159 — 00000045 | .00000014 .00000001 

12.0 | -00000000 | .00000099 — .00000003 | .0000001 1 — .00000004 

12.5 | .00003231 | —.00000007 | 00000012 | .00000002 — .00000002 

13.0 | .00000000 — .00000034 | .00000007 .00000001 — .00000001 

13.5 | —.00001492 — .00000014 | .00000000 .00000001 

14.0 | .00000000 00000004 | —.00000002 

14.5 | .00000690 | 00000007 | —.00000001 

15.0 .00000000 .00000002 

15.5 — 00000318 | —.00000001 | 

16.0 .00000000 | —.00000001 

16.5 | .00000147 | 

17.0 | .00000000 | 

17.5 | —.00000068 | 

18.0 | .00000000 | 

18.5 | .00000031 | | 

19.0 | .00000000 | | 

19.5 | —.00000014 | | 

20.0 | 00000000 | 

20.5 | .00000007 | | 

21.0 | .00000000 | | 

21.5 — .00000003 | 

22.0 00000000 | | 

72.3 | .00000001 | 

23.0 | .00000000 | | 

23.5 - 


| 


00000001 
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.00000001 | 


























L,(x, 1/2, ©) 
e=.5 e=.6 e=x.7 | e=.8 | e=.9 e=1.0 
.59702260 | .58765637 | .58031608 .57438799 | .56948898 .56536580 
.47675954 47077398 | .46599416 .46207717 .45880202 .45601892 
23077657 | .23140004 | . 23168090 | .23177314 .23175789 .23168050 
.04557847 .05027848 -05380648 | -05654957 | .05874137 .06053136 
.02276409 | —.01820178 .01459584 | 01167397 | .00925829 .00722771 
.02384227 | — .02167108 .01979233 | .01816756 | .01675611 .01552233 
- 00896161 | — .00898986 | .00881784 | -00855221 .00824659 .00792883 
.00044979 — .00038991 .00093726 .00129962 .00154092 .00170083 
.00238588 | .00180229 .00135734 -00101562 .00075048 .00054256 
.00127570 | .00114308 -00100305 | .00087278 | .00075725 .00065681 
.00019600 .00027865 | .00030985 | .00031471 .00030617 .00029109 
.00019071 | —.00009654 .00003598 | .00000283 | .00002756 .00004311 
.00015995 | — .00012137 | .00008970 .00006510 .00004639 .00003223 
.00004696 | —.00004885 | .00004474 | .00003886 | .00003288 .00002744 
.00000987 | —.00000179 | 00000738 | .00000973 .00001039 .00001018 
.00001687 .00001000 | .00000536 00000238 | .00000050 .00000064 
.00000771 .00000641 | .00000486 .00000350 .00000244 .00000165 
00000044 | .00000148 .00000168 .00000156 .00000134 .00000110 
00000156 | —.00000057 | .00000004 .00000021 .00000031 .00000034 
.00000101 | —.00000067 .00000039 -00000020 | .00000009 .00000002 
.00000023 | — .00000027 | .00000022 | .00000017 | .00000011 | .00000008 
.00000011 | —.00000001 | .00000005 | -00000006 | .00000005 .00000004 
.00000012 .00000006 | .00000002 .00000000 | .00000001 | .00000001 
.00000004 | -00000004 .00000002 .00000001 .00000001 
.00000000 .00000001 .00000001 | .00000001 
.00000001 | 
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Li! (x, 1/2, ©) 

x e=.0 e=,1 | e=,.2 e=.3 e=.4 
0.0 —3.47753764 —1.50781449 | —1.27083552 —1.16381926 —1.10100908 
0.5 —1.03983382 — .69689346 — .61542693 — .57326418 — .54670549 
1.0 2.15613767 .54167607 | .40225158 .34798919 .31925128 
1 .50655095 .77131824 -63355043 .56889199 .53067521 
2.0 — .58680458 .31875073 | .30878327 . 29072626 .27601758 
z5 — .68097024 — .03482567 .02460402 -04246822 .04943501 
3.0 . 24590776 — .12647279 — .07651577 — .05154399 — .03726679 
.31311773 — .06455381 — .05654901 — .04581150 — .03775297 
4.0 — .11165611 .01678358 — .00530735 — .01047426 — .01153447 
4.5 — .14452585 .03142765 .01425664 .00665986 .00297134 
5.0 .05142591 .00673809 .00811324 -00596871 .00423710 
aS .06675341 — .00655085 .00060896 .00183704 .00187944 
6.0 — .02374375 — .00477107 — .00194355 — .00054832 .00001970 
6.5 — .03083551 — .00059653 — .00138896 — .00087620 — .00050023 
7.0 .01096729 .00133431 — .00011460 — .00030729 — .00026914 
ta .01424417 .00097571 | -00035809 .00005960 — .00003417 
8.0 — .00506618 — .00005726 | .00019854 .00010609 .00004375 
8.5 — .00657998 — .00035828 | .00001209 .00004959 .00003703 
9.0 .00234028 — .00012112 | — .00004884 — .00000116 .00000954 
9.5 .00303957 .00004878 | — .00003391 — .00001505 — .00000404 

10.0 — .00108107 .00005884 | — .00000242 — .00000732 — .00000433 
10.5 — .00140411 .00001802 | .00000898 — .00000026 — .00000164 
11.0 .00049939 — .00001229 .00000485 .00000168 .00000013 
5 .00064862 — .00001316 .00000023 .00000113 .00000053 
12.0 — .00023069 — .00000110 — .00000123 .00000016 -00000025 
12.5 — .00029962 .00000374 — .00000083 — .00000023 .00000001 
13.0 .00010657 .00000182 — .00000005 — .00000015 — .00000005 
13.5 .00013841 — .00000018 .00000022 — .00000003 — .00000004 
14.0 — .00004923 — .00000067 .00000012 .00000002 — .00000001 
14.5 — .00006394 — .00000033 .00000000 -00000002 .00000001 
15.0 .00002274 -00000009 — .00000003 .00000001 

15.5 .00002954 .00000016 — .00000002 

16.0 — .00001050 .00000003 

16.5 — .00001364 — .00000003 

17.0 .00000485 — .00000002 

i .00000630 .00000000 

18.0 — .00000224 .00000001 

18.5 — .00000291 

19.0 .00000104 

19.5 .00000134 

20.0 — .00000048 

20.5 — .00000062 

21.0 .00000022 

2.9 .00000029 

22.0 — .00000010 

0 — .00000013 

23.0 .00000005 

23.3 -00000006 

24.0 — .00000002 

24.5 — .00000003 

25.0 .00000001 

20.9 .00000001 

26.0 — .00000001 
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Li (x, 1/2, &) 























x e=.5 e=.6 e=.7 e=.8 e=.9 e=1.0 
0.0 |—1.05919265 |—1.02916420 |—1.00647330 — .98868331 | —.97433986 | —.96251767 
0.5 | —.52821280 | —.51450880 | —.50390754 | —.49544281 | —.48851722 | —.48273978 
1.0 30155959 . 28962736 | 28106624 27464186 . 26965348 . 26567454 
1.5 50527188 48711044 | .47346009 46281678 45428116 44728133 
2.0 26453429 . 25546310 | . 24815810 24216442 23716440 . 23293281 
2.5 | .05240385 05363758 |  .05404202 .05402794 .05379802 05345832 
3.0 — .02823788 | —.02210889 | —.01772635 | —.01446547 | —.01196172 | —.00998973 
3.5 — 03182937 | —.02737565 | —.02393662 | —.02121579 | —.01901742 | —.01720888 
4.0 — .01135624 | —.01078633 | — .01011492 | —.00944820 | —.00882438 | —.00825481 
4.5 .00101034 | —.00009918 | —.00075259 | — .00114593 | — .00138383 | —.00152531 
5.0 | 00303129 00219629 |  .00160744 .00118299 .00087064 00063651 
5.5 | 00166647 00142136 .00120000 00101295 00085821 .00073079 
6.0 | .00024844 00033324 |  .00035457 .00034759 00032845 .00030479 
6.5 — .00027461 | —.00014145 | —.00006202 | — .00001410 00001495 .00003245 
7.0 — .00020347 | —.00014764 | —.00010563 | —.00007505 | —.00005291 | —.00003686 
7.5 — .00005788 | —.00005841 | —.00005187 | — .00004389 | — .00003637 | — .00002987 
8.0 | 00001265 | —.00000184 | — .00000816 | — .00001051 | —.00001095 | — .00001051 
8.5 00002301 00001329 | .00000717 .00000344 .00000121 | —.00000014 
9.0 | .00000979 .00000777 | .00000567 .00000399 00000274 00000184 
9.5 | .00000025 .00000163 |  .00000187 00000172 00000145 .000001 17 
10.0 — .00000199 | —.00000071 | — .00000008 00000021 .00000032 .00000034 
10.5 — .00000133 | —.00000085 | —.00000049 | — .00000026 | —.00000012 | — .00000004 
11.0 | —.00000029 | —.00000033 | — .00000026 | —.00000019 | —.00000013 | — .00000008 
11.5 |  .00000016 00000000 | — .00000005 | — .00000006 | —.00000005 | — .00000004 
12.0 | 00000015 .00000007 | .00000002 00000000 | —.00000001 | — .00000001 
12.5 | 00000005 .00000004 |  .00000003 .00000001 .00000001 
13.0 |  .00000000 .00000001 | .00000001 00000001 
13.5 — .00000002 
14.0 — .00000001 
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—NOTES— 
SOME APPLICATIONS OF THE REPEATED INTEGRALS 
OF THE ERROR FUNCTION* 
By J. C. JAEGER (University of Tasmania) 


1. Introductory. The repeated integrals of the error function 


i" erfc x = f i" erfc édé, c= 1, 2,«+- (1) 


where 


9 Ce) 
i° erfc x = erfc x = f e * dé, (2) 
VTH x2 


have been studied by Hartree,! who tabulates them for »=1 and n=2, and shows 
that they satisfy the recurrence relation 


2n i* erfc x = i** erfe x — 2xi"" erfc x. (3) 
He also shows that 
L{ (44)"/?i" erfc (Zat-!/2)} = s—I-n/2e-eve, (4) 
where 2=0,1,2,---,a20,andZ fv} is written for the Laplace transform of a func- 
tion v(t) of ¢, that is ‘ 
L{v} = f e~*'y(t)dt. (5) 
0 


The functions (1) arise naturally in the theory of conduction of heat in the semi- 
infinite solid (or the sphere or slab) with prescribed surface temperature or flow of 
heat, since the Laplace transforms of the solutions of many such problems involve 
the functions on the right hand side of (4). 

The objects of this note are, firstly, to indicate an extension of (4) which applies 
in the same way to problems with heat transfer at the surface, and, secondly, to give 
solutions in terms of the functions (1) of a number of problems of practical interest 
which involve heat generation in the solid. 

2. Problems involving heat transfer at a surface at a rate proportional to its 
temperature difference from its surroundings. The required extension of (4) is 
that, if ” is a positive integer, and a, h, and x are positive, 


n—1 —qz 
thas h)-* | eon erfc (H + X) — >> (— 2H)’i" erfc x | = pein (6) 
\ ane g**(q -t- h) 
where 
X=2/2Vat, H=hJat, q=Vvs/a. (7) 


* Received Sept. 19, 1945. 
1D. R. Hartree, Some properties and applications of the repeated integrals of the error function, Proc. 
Manchester Lit. and Phil. Soc. 80, 85 (1935). 
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To derive this result we notice that 


eu (—)*"e~9 eu n—1 h r 
Soa + oen(-). (8) 
ghi(qt+h) heqgqt+h) hq? 0 q 


In the terms of the series we use (4); the result for the first term of the right-hand 
side of (8) is given in most tables of Laplace transforms; and (6) follows immediately. 

Typical examples in which (6) arises are the following: 

(i) The semi-infinite solid x>0. Zero initial temperature. The solid heated at x =0 
for t>0 by heat transfer from a medium at at"’?, n=0,1,---. 

The temperature v in the solid has to satisfy 








i cp ae Ae (9) 


with the boundary condition 


Ov 
— = — h(at*/? — »), z=0, #>0. 
Ox 
Also v has to be bounded as x ©. The Laplace transform of the solution is found 


to be 


haY (1 + n/2)e~97 
sitnl2(g + h) 
using the notation (7). Therefore, from (6) 
(—)**!aP(1 + n/2) ( ® 
v= SI ih Ss ad erfc (H + X) — >> (— 2A)ri’ erfc x} ; 
r=0 


h*a”!? 


L{v} = 





If heat transfer takes place from a medium whose temperature is 
N 
> ant"! (10) 
n=0 


the solution follows at once. For the problem of the semi-infinite solid whose surface 
temperature is given by (10) the result is obtained in the same way by using (4). 
An empirical relation of the type (10) is often useful for representing observed surface 
temperatures, the term in t'/? being particularly valuable since it corresponds to con- 
stant flux of heat; for example the fall in temperature of the Earth’s surface after 
sunset on a cloudless evening is approximately proportional to ¢t'/? and may be repre- 
sented very well by two or three terms of (10). 

(ii) The semi-infinite solid x>0, of conductivity K and diffusivity a. At x=0 the 
solid is in contact with mass M per unit area of well stirred fluid of specific heat c’, 
whose temperature is equal to the surface temperature of the solid, i.e. to limz.40v. The int- 
tial temperatures of the solid and fluid are zero. Heat is supplied to the fluid at constant 
rate Q per unit mass per unit time for t>0. 

Here (9) has to be solved with boundary condition at x =0 


Ov ov 
Mc' —-— K— = QM, 
ot Ox 
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The Laplace transform of the solution is 

Qe 
c'g°(q + h) 
where now h=K/Mc’a. Thus from (6) we have 


evo 
v= = wy { e2H 
h?ac 





L{v} = 


X+X* erfc (H + X) — erfe X + 2Hi! erfc X}. 


(iii) The semi-infinite solid x >0. Zero initial temperature. Heat is produced for t>0 
in the solid at the rate Qi"/?, n= —1,0,1,---, per unit time per unit volume. There is 
heat transfer at x =0 into a medium at zero temperature. 

Here we have to solve the differential equation 

079 1 dv Q 
——-— — = — —#"% x > 0, #>0, 
Ox? a dt K 


with boundary condition 


Ov ' 
—— hv=0, g= 0, 2>°0. 
Ox 
Here 
iano jan a... 2 
K gttn/2 a2tnl2gnt4(g a h)§ 
and 
atitn/2 IT (1 + n/2) . at? 
,= _ ee Ht Bid ds — | sono erfc (H + X) — 2% (— 2H)’i’ erfc x |. 
K(1 + n/2) Ka*!2(— h)**? ser 


3. Cases of generation of heat in a solid. The solutions of a number of problems 
of practical importance in which heat is generated for ¢>0 in a solid at the rate 


N 


>> ant”? (11) 


n=—l 


per unit time per unit volume can be expressed in terms of the functions (1). An ex- 
pression of type (11) may be useful for representing an experimentally observed rate 
of generation of heat; the term in ¢~/? is of value when the initial rate of heat produc- 
tion is high, as in the hydrating of cement.? 

A problem involving (11) and radiation at the surface has already been given in 
§2(iii), here we give the solutions of some cases in which the surface temperature is 
zero. 

(iv) The region x>0. x=0 kept at zero for t>0. Heat production at the rate t"!* in 
0 <x <a, and zero in x>a. Zero initial temperature. 


2 In this case two or three terms of (11) give quite a good representation of the observed heat of hy- 
dration over the first few days. An expression in terms of negative exponentials is more usual, but does 
not lead to solutions in terms of tabulated functions for the problems given here, except (vi), and in that 
case it does not give a solution which is useful for small values of the time. 
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This is the fundamental practical problem of temperatures in hydrating concrete: 
a slab of concrete is poured on the surface of the semi-infinite solid, we assume here 
that the thermal constants of the concrete and the solid are the same. The solution is 


afitels J : a-—2z , a+x 
y= ————-_ 41 — ['(2 + n/2)2"*'| i? erfc ———— — i**? erfc ———— 
K(1 + n/2) \ 2(at)!/2 2(at)"/2 


+ 2i"*? erfc sau |} 


2(at)!/? 
if 0<x <a, and for x><a it is 
al'(1 + n/2)(4t)'*"/? ¢. z—a ; x+a 
,= ——_—— {it? erfe ———— + i**? erfc ——— 
2K 2(at)*/? 2(at)/? 


x 
— 2i"*? erfc sant . 
2(at) 1/2 
(v) The problem* of (iv) except that heat is produced only in the region a<x <b. 


The temperature gradient at the surface is 


ql? a b 
— I(1 + n/2)(4t a Lh erfc ————  — i'** erfc san | 
K ) 2(at)!/? 2(at)*/? 


(vi) The slab 0<x<l. The surfaces x=0 and x=l kept at zero temperature. Heat 


generation at the rate t"!*. Zero initial temperature. 


afitn!2 J “ ml + x 
v= 1 — 1(2 + n/2)2"*? — -| it? erfc - — 
K(1 + n/2) \ 2 2(at)'/? 
' (m+ 1)l—x 
+ i*+? erf¢c ————_ \ ‘ 
2(at)*/? 


(vii) The infinite region r=0. Zero initial temperature. Heat production at the rate 


t"/2 in the sphere 0 Sr <a, zero elsewhere. 











afi+n/2 aal(1 + n/2)(40)'*"/? a-?r ; a+r 
v= — = [in erfc — — i*+? erfc ———_ 
K(1 + n/2) 2Kr 2(at) 1/2 2(at)1/? 
aF/?T(1 + n 2) (4t) (nt8/2 A a= A a+r 
- - ——| i*** eric ———— — i*** erfc — —|,0Sr<a. 
2Kr 2(at)'!? 2(at)*'/? 
aal'(1 + n/2)(40)'*"?7 | a+r ; r— <a 
v= _ ——j} iv? eric —— > “ete ——— 
2Kr 2(at)1/? 2(at)*/? 
aF/2P(1 + n/2)(4t) rt ?? a+r ; r—a@ 
——_——_——— | erfc ———— — i*** erfc ——— |], r>a. 
2Kr 2(at) 1/2 2(axt) 1/2 
3 Van Ostrand, On the flow of heat from a rock stratum in which heat is being generated, J. Wash. Acad. 
Sci. 22, 529 (1932), considers the case of a thin layer. 


4 This problem is of interest in connection with development of heat in wheat stacks. 
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The development of mathematics. By E. T. Bell. Second edition. McGraw-Hill Book 
Company, Inc., New York and London, 1945. xiii+637 pp. $5.00. 


The second edition of this well known book contains about fifty pages of new material. Listed in the 
order in which they occur in the book, the major additions deal with Moslem algebra, the development 
of symbolism, lattice theory, cubic surfaces, Levi-Civita’s parallel displacement, definition of lengths, 
areas and volumes, the transition from intuitive to unintuitive thinking in modern mathematics, the de- 
velopment of mathematics in times of war, algebra of relations and consistency proofs. The index, too, 
has been considerably enlarged. Shorter additional paragraphs are concerned with applied mathematics 
in World War II, Egyptian algebra, number mysticism, Greek methods of computation, the Egyptian 
method of constructing a right angle, mathematical realism, the Greek treatment of loci, number theory, 
geometrical constructions, projective differential geometry, theory of quantics, intuitionism, Whitehead 
and Russell’s Principia Mathematica, the Burali-Forti paradox, quantum theory, Skolem’s theorem, 


Laplace’s work on probability, R. A. Fisher’s work on statistics, and matrix algebra in modern statistics, 
S. PRAGER 


Tables of associated Legendre functions. Prepared by the Mathematical Tables Proj- 
ect conducted under the sponsorship of the National Bureau of Standards. Offi- 
cial Sponsor: Lyman J. Briggs. Project Director: Arnold N. Lcwan. Columbia 
University Press, New York, 1945. xxv+303 pp. $5.00. 


The main body of the present volume contains tables for P™(x), dP™(x)/dx, 4-"P™ (ix), i-"dP (ix) /dx, 
(—1)™Qm™(x), (—1)™*4dQ™(x)/dx, —1"t2™+1Om (ix), ant2™—1dQm (ix) /dx, Prysj2(x), dPry1/2(x)/dx, (—1)™Qr 12x) 
and (—1)"*dQ%.1/2(x)/dx for integral values of m (ranging from 0 to 4) and n (ranging from 0 or 1 to 
10 for the functions of integral degree and from —1 to 4 for the functions of half-integral degree). For the 
greater part, these functions are tabulated to six significant figures for values of the argument which in- 
crease from 0 or 1 to 10 by steps of 0.1. The functions P™(6) and dP™(6)/dx are also tabulated for integral 
values of m (ranging from 0 or 1 to 4) and n (ranging from 1 to 10), the values of the argument increasing 


from 0° to 90° by steps of 1°. Auxiliary tables facilitate interpolation. 
W. PRAGER 
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